1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <ctype.h>
5 #include <math.h>
6
7 #include "gradsdes.h"
8
9 static char pout[1024];
10 FILE *descr; /* File descriptor pointer */
11 int cal365 = -1;
12 int fullyear = -999;
13
14 void
dsets_init(dsets_t * pfi)15 dsets_init(dsets_t *pfi)
16 {
17 pfi->name[0] = 0;
18 pfi->dnam[0] = 0;
19 pfi->title[0] = 0;
20 pfi->bswap = 0;
21 pfi->fhdr = 0;
22 pfi->xyhdr = 0;
23 pfi->seqflg = 0;
24 pfi->yrflg = 0;
25 pfi->zrflg = 0;
26 pfi->flt64 = 0;
27 pfi->tmplat = 0;
28 pfi->pa2mb = 0;
29 pfi->calendar = 0;
30 pfi->type = 1; /* Assume grid unless told otherwise */
31 pfi->idxflg = 0; /* Assume binary */
32 pfi->ncflg = 0; /* Assume not netcdf */
33
34 pfi->undef = -9.99E33;
35
36 pfi->pvar1 = NULL;
37 pfi->ens1 = NULL;
38
39 pfi->pchsub1 = NULL;
40
41 for (int i = 0; i < 5; ++i) pfi->dnum[i] = 0;
42 }
43
44 /* Byte swap requested number of 4 byte elements */
45
46 void
gabswp(void * r,gaint cnt)47 gabswp(void *r, gaint cnt)
48 {
49 gaint i;
50 char *ch1, *ch2, *ch3, *ch4, cc1, cc2;
51
52 ch1 = (char *) r;
53 ch2 = ch1 + 1;
54 ch3 = ch2 + 1;
55 ch4 = ch3 + 1;
56 for (i = 0; i < cnt; i++)
57 {
58 cc1 = *ch1;
59 cc2 = *ch2;
60 *ch1 = *ch4;
61 *ch2 = *ch3;
62 *ch3 = cc2;
63 *ch4 = cc1;
64 ch1 += 4;
65 ch2 += 4;
66 ch3 += 4;
67 ch4 += 4;
68 }
69 }
70
71 /* Byte swap requested number of 2 byte elements */
72
73 void
gabswp2(void * r,gaint cnt)74 gabswp2(void *r, gaint cnt)
75 {
76 gaint i;
77 char *ch1, *ch2, cc1;
78
79 ch1 = (char *) r;
80 ch2 = ch1 + 1;
81 for (i = 0; i < cnt; i++)
82 {
83 cc1 = *ch1;
84 *ch1 = *ch2;
85 *ch2 = cc1;
86 ch1 += 2;
87 ch2 += 2;
88 }
89 }
90
91 /*mf version
92 convert all upper case alphabetic characters to lower case.
93 The GrADS system is case insensitive, and assumes lower case
94 internally in most cases. Does not turn to lower case if in "'s
95 */
96 void
lowcas(char * ch)97 lowcas(char *ch)
98 {
99 int i;
100 int qflag = 0;
101
102 while (*ch != '\0' && *ch != '\n')
103 {
104 i = *ch;
105 if (*ch == '\"' && qflag == 0)
106 {
107 qflag = 1;
108 }
109 else if (*ch == '\"' && qflag == 1)
110 {
111 qflag = 0;
112 }
113 if (i > 64 && i < 91 && qflag == 0)
114 {
115 i += 32;
116 *ch = i;
117 }
118 else if (i == 95)
119 {
120 *ch = i;
121 }
122 ch++;
123 }
124 }
125
126 /* Date/Time manipulation routines. Note that these routines
127 are not particularly efficient, thus Date/Time conversions
128 should be kept to a minimum. */
129
130 static gaint mosiz[13] = { 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
131 static gaint momn[13] = { 0, 44640, 40320, 44640, 43200, 44640, 43200, 44640, 44640, 43200, 44640, 43200, 44640 };
132 static gaint mnacum[13] = { 0, 0, 44640, 84960, 129600, 172800, 217440, 260640, 305280, 349920, 393120, 437760, 480960 };
133 static gaint mnacul[13] = { 0, 0, 44640, 86400, 131040, 174240, 218880, 262080, 306720, 351360, 394560, 439200, 482400 };
134
135 /* Test for leap year. Rules are:
136
137 Divisible by 4, it is a leap year, unless....
138 Divisible by 100, it is not a leap year, unless...
139 Divisible by 400, it is a leap year. */
140
141 gaint
qleap(gaint year)142 qleap(gaint year)
143 {
144 gaint i, y;
145
146 /*mf - disable if 365 day calendar mf*/
147
148 if (/*mfcmn.*/ cal365 == 1) return (0);
149
150 y = year;
151
152 i = y / 4;
153 i = (i * 4) - y;
154 if (i != 0) return 0;
155
156 i = y / 100;
157 i = (i * 100) - y;
158 if (i != 0) return (1);
159
160 i = y / 400;
161 i = (i * 400) - y;
162 if (i != 0) return 0;
163
164 return (1);
165 }
166
167 /* Add an offset to a time. Output to dto. */
168
169 void
timadd(struct dt * dtim,struct dt * dto)170 timadd(struct dt *dtim, struct dt *dto)
171 {
172 gaint i;
173 gaint cont;
174
175 /* First add months and years. Normalize as needed. */
176 dto->mo += dtim->mo;
177 dto->yr += dtim->yr;
178
179 while (dto->mo > 12)
180 {
181 dto->mo -= 12;
182 dto->yr++;
183 }
184
185 /* Add minutes, hours, and days directly. Then normalize
186 to days, then normalize extra days to months/years. */
187
188 dto->mn += dtim->mn;
189 dto->hr += dtim->hr;
190 dto->dy += dtim->dy;
191
192 if (dto->mn > 59)
193 {
194 i = dto->mn / 60;
195 dto->hr += i;
196 dto->mn = dto->mn - (i * 60);
197 }
198 if (dto->hr > 23)
199 {
200 i = dto->hr / 24;
201 dto->dy += i;
202 dto->hr = dto->hr - (i * 24);
203 }
204
205 cont = 1;
206 while (dto->dy > mosiz[dto->mo] && cont)
207 {
208 if (dto->mo == 2 && qleap(dto->yr))
209 {
210 if (dto->dy == 29)
211 cont = 0;
212 else
213 {
214 dto->dy -= 29;
215 dto->mo++;
216 }
217 }
218 else
219 {
220 dto->dy -= mosiz[dto->mo];
221 dto->mo++;
222 }
223 while (dto->mo > 12)
224 {
225 dto->mo -= 12;
226 dto->yr++;
227 }
228 }
229 }
230
231 /* Subtract an offset from a time. Subtract minutes/hours/days
232 first so that we will exactly reverse the operation of timadd */
233
234 void
timsub(struct dt * dtim,struct dt * dto)235 timsub(struct dt *dtim, struct dt *dto)
236 {
237 gaint s1, s2;
238
239 /* Subtract minutes, hour, and days directly. Then normalize
240 to days, then normalize deficient days from months/years. */
241
242 dto->mn = dtim->mn - dto->mn;
243 dto->hr = dtim->hr - dto->hr;
244 dto->dy = dtim->dy - dto->dy;
245 s1 = dto->mo;
246 s2 = dto->yr;
247 dto->mo = dtim->mo;
248 dto->yr = dtim->yr;
249
250 while (dto->mn < 0)
251 {
252 dto->mn += 60;
253 dto->hr--;
254 }
255 while (dto->hr < 0)
256 {
257 dto->hr += 24;
258 dto->dy--;
259 }
260
261 while (dto->dy < 1)
262 {
263 dto->mo--;
264 if (dto->mo < 1)
265 {
266 dto->mo = 12;
267 dto->yr--;
268 }
269 if (dto->mo == 2 && qleap(dto->yr))
270 dto->dy += 29;
271 else
272 dto->dy += mosiz[dto->mo];
273 }
274
275 /* Now subtract months and years. Normalize as needed. */
276
277 dto->mo = dto->mo - s1;
278 dto->yr = dto->yr - s2;
279
280 while (dto->mo < 1)
281 {
282 dto->mo += 12;
283 dto->yr--;
284 }
285
286 /* Adjust for leaps */
287
288 if (dto->mo == 2 && dto->dy == 29 && !qleap(dto->yr))
289 {
290 dto->mo = 3;
291 dto->dy = 1;
292 }
293 }
294
295 /* Convert from a t grid coordinate to an absolute time. */
296
297 void
gr2t(gadouble * vals,gadouble gr,struct dt * dtim)298 gr2t(gadouble *vals, gadouble gr, struct dt *dtim)
299 {
300 struct dt stim;
301 gadouble *moincr, *mnincr;
302 gadouble v;
303
304 /* Get constants associated with this conversion */
305 stim.yr = (gaint)(*vals + 0.1);
306 stim.mo = (gaint)(*(vals + 1) + 0.1);
307 stim.dy = (gaint)(*(vals + 2) + 0.1);
308 stim.hr = (gaint)(*(vals + 3) + 0.1);
309 stim.mn = (gaint)(*(vals + 4) + 0.1);
310 moincr = vals + 5;
311 mnincr = vals + 6;
312
313 /* Initialize output time */
314 dtim->yr = 0;
315 dtim->mo = 0;
316 dtim->dy = 0;
317 dtim->hr = 0;
318 dtim->mn = 0;
319
320 /* Do conversion if increment is in minutes. */
321 if (*mnincr > 0.1)
322 {
323 v = *mnincr * (gr - 1.0);
324 if (v > 0.0)
325 v = v + 0.5; /* round */
326 else
327 v = v - 0.5;
328 dtim->mn = (gaint) v;
329 if (dtim->mn < 0)
330 {
331 dtim->mn = -1 * dtim->mn;
332 timsub(&stim, dtim);
333 }
334 else
335 {
336 timadd(&stim, dtim);
337 }
338 return;
339
340 /* Do conversion if increment is in months. Same as for minutes,
341 except special handling is required for partial months.
342 JMA There is a bug here, and some precision decisions that need
343 attention */
344 }
345 else
346 {
347 v = *moincr * (gr - 1.0);
348 if (v < 0.0)
349 dtim->mo = (gaint)(v - 0.9999); /* round (sort of) */
350 else
351 dtim->mo = (gaint)(v + 0.0001);
352 v = v - (gadouble) dtim->mo; /* Get fractional month */
353 if (dtim->mo < 0)
354 {
355 dtim->mo = -1 * dtim->mo;
356 timsub(&stim, dtim);
357 }
358 else
359 timadd(&stim, dtim);
360 if (v < 0.0001) return; /* if fraction small, return */
361
362 if (dtim->mo == 2 && qleap(dtim->yr))
363 {
364 v = v * 41760.0;
365 }
366 else
367 {
368 v = v * (gadouble) momn[dtim->mo];
369 }
370 stim = *dtim;
371 dtim->yr = 0;
372 dtim->mo = 0;
373 dtim->dy = 0;
374 dtim->hr = 0;
375 dtim->mn = (gaint)(v + 0.5);
376 timadd(&stim, dtim);
377 return;
378 }
379 }
380
381 /* Calculate the difference between two times and return the
382 difference in minutes. The calculation is time2 - time1, so
383 if time2 is earlier than time1, the result is negative. */
384
385 gaint
timdif(struct dt * dtim1,struct dt * dtim2)386 timdif(struct dt *dtim1, struct dt *dtim2)
387 {
388 gaint min1, min2, yr;
389 struct dt *temp;
390 gaint swap, mo1, mo2;
391
392 swap = 0;
393 if (dtim1->yr > dtim2->yr)
394 {
395 temp = dtim1;
396 dtim1 = dtim2;
397 dtim2 = temp;
398 swap = 1;
399 }
400
401 min1 = 0;
402 min2 = 0;
403
404 yr = dtim1->yr;
405 while (yr < dtim2->yr)
406 {
407 if (qleap(yr))
408 min2 += 527040L;
409 else
410 min2 += 525600L;
411 yr++;
412 }
413
414 mo1 = dtim1->mo;
415 mo2 = dtim2->mo;
416 if (qleap(dtim1->yr))
417 {
418 min1 = min1 + mnacul[mo1] + (dtim1->dy * 1440L) + (dtim1->hr * 60L) + dtim1->mn;
419 }
420 else
421 {
422 min1 = min1 + mnacum[mo1] + (dtim1->dy * 1440L) + (dtim1->hr * 60L) + dtim1->mn;
423 }
424 if (qleap(dtim2->yr))
425 {
426 min2 = min2 + mnacul[mo2] + (dtim2->dy * 1440L) + (dtim2->hr * 60L) + dtim2->mn;
427 }
428 else
429 {
430 min2 = min2 + mnacum[mo2] + (dtim2->dy * 1440L) + (dtim2->hr * 60L) + dtim2->mn;
431 }
432 if (swap)
433 return (min1 - min2);
434 else
435 return (min2 - min1);
436 }
437
438 static const char *mons[12] = { "jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec" };
439
440 /* Parse an absolute date/time value. Format is:
441
442 12:00z 1jan 1989 (jan,feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec)
443
444 Must have Z or Month abbrev, or value is invalid. 'def' contains
445 higher order missing values (usually from tmin in pst). Lower order
446 values are defaulted to be: dy = 1, hr = 0, mn = 0. */
447
448 char *
adtprs(char * ch,struct dt * def,struct dt * dtim)449 adtprs(char *ch, struct dt *def, struct dt *dtim)
450 {
451 gaint val, flag, i;
452 char monam[5];
453
454 dtim->mn = 0;
455 dtim->hr = 0;
456 dtim->dy = 1;
457
458 if (*ch >= '0' && *ch <= '9')
459 {
460 flag = 0;
461 ch = intprs(ch, &val);
462 if (*ch == ':' || tolower(*ch) == 'z')
463 {
464 if (val > 23)
465 {
466 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
467 sprintf(pout, " Hour = %i -- greater than 23\n", val);
468 gaprnt(0, pout);
469 return (NULL);
470 }
471 dtim->hr = val;
472 if (*ch == ':')
473 {
474 ch++;
475 if (*ch >= '0' && *ch <= '9')
476 {
477 ch = intprs(ch, &val);
478 if (val > 59)
479 {
480 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
481 sprintf(pout, " Minute = %i -- greater than 59\n", val);
482 gaprnt(0, pout);
483 return (NULL);
484 }
485 if (tolower(*ch) != 'z')
486 {
487 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
488 gaprnt(0, " 'z' delimiter is missing \n");
489 return (NULL);
490 }
491 dtim->mn = val;
492 ch++;
493 if (*ch >= '0' && *ch <= '9')
494 ch = intprs(ch, &val);
495 else
496 val = def->dy;
497 }
498 else
499 {
500 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
501 gaprnt(0, " Missing minute value \n");
502 return (NULL);
503 }
504 }
505 else
506 {
507 ch++;
508 if (*ch >= '0' && *ch <= '9')
509 ch = intprs(ch, &val);
510 else
511 val = def->dy;
512 }
513 }
514 else
515 flag = 2;
516 dtim->dy = val;
517 }
518 else
519 flag = 1;
520
521 monam[0] = tolower(*ch);
522 monam[1] = tolower(*(ch + 1));
523 monam[2] = tolower(*(ch + 2));
524 monam[3] = '\0';
525
526 i = 0;
527 while (i < 12 && !cmpwrd(monam, mons[i])) i++;
528 i++;
529
530 if (i == 13)
531 {
532 if (flag == 1)
533 {
534 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
535 gaprnt(0, " Expected month abbreviation, none found\n");
536 return (NULL);
537 }
538 if (flag == 2)
539 {
540 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
541 gaprnt(0, " Missing month abbreviation or 'z' delimiter\n");
542 return (NULL);
543 }
544 dtim->mo = def->mo;
545 dtim->yr = def->yr;
546 }
547 else
548 {
549 dtim->mo = i;
550 ch += 3;
551 /* parse year */
552 if (*ch >= '0' && *ch <= '9')
553 {
554 /* use fullyear only if year 1 = 0001*/
555 if (*(ch + 2) >= '0' && *(ch + 2) <= '9')
556 {
557 /*mfcmn.*/ fullyear = 1;
558 }
559 else
560 {
561 /*mfcmn.*/ fullyear = 0;
562 }
563 ch = intprs(ch, &val);
564 }
565 else
566 {
567 val = def->yr;
568 }
569
570 /* turn off setting of < 100 years to 1900 or 2000 */
571 if (/*mfcmn.*/ fullyear == 0)
572 {
573 if (val < 50)
574 val += 2000;
575 else if (val < 100)
576 val += 1900;
577 }
578 dtim->yr = val;
579 }
580
581 i = mosiz[dtim->mo];
582 if (dtim->mo == 2 && qleap(dtim->yr)) i = 29;
583 if (dtim->dy > i)
584 {
585 gaprnt(0, "Syntax Error: Invalid Date/Time value.\n");
586 sprintf(pout, " Day = %i -- greater than %i \n", dtim->dy, i);
587 gaprnt(0, pout);
588 return (NULL);
589 }
590 return (ch);
591 }
592
593 /* Parse a relative date/time (offset). Format is:
594
595 nn (yr/mo/dy/hr/mn)
596
597 Examples: 5mo
598 1dy12hr
599 etc.
600
601 Missing values are filled in with 0s. */
602
603 char *
rdtprs(char * ch,struct dt * dtim)604 rdtprs(char *ch, struct dt *dtim)
605 {
606 gaint flag, val;
607 char id[3];
608
609 dtim->yr = 0;
610 dtim->mo = 0;
611 dtim->dy = 0;
612 dtim->hr = 0;
613 dtim->mn = 0;
614
615 flag = 1;
616
617 while (*ch >= '0' && *ch <= '9')
618 {
619 flag = 0;
620 ch = intprs(ch, &val);
621 id[0] = *ch;
622 id[1] = *(ch + 1);
623 id[2] = '\0';
624 if (cmpwrd("yr", id))
625 dtim->yr = val;
626 else if (cmpwrd("mo", id))
627 dtim->mo = val;
628 else if (cmpwrd("dy", id))
629 dtim->dy = val;
630 else if (cmpwrd("hr", id))
631 dtim->hr = val;
632 else if (cmpwrd("mn", id))
633 dtim->mn = val;
634 else
635 {
636 gaprnt(0, "Syntax Error: Invalid Date/Time offset.\n");
637 sprintf(pout, " Expecting yr/mo/dy/hr/mn, found %s\n", id);
638 gaprnt(0, pout);
639 return (NULL);
640 }
641 ch += 2;
642 }
643 if (flag)
644 {
645 gaprnt(0, "Syntax Error: Invalid Date/Time offset.\n");
646 gaprnt(0, " No offset value given\n");
647 return (NULL);
648 }
649 return (ch);
650 }
651
652 /* Compares two strings. A match occurs if the leading
653 blank-delimited words in the two strings match. CR and NULL also
654 serve as delimiters. */
655
656 gaint
cmpwrd(const char * ch1,const char * ch2)657 cmpwrd(const char *ch1, const char *ch2)
658 {
659
660 while (*ch1 == ' ' || *ch1 == '\t') ch1++; /* Advance past leading blanks. */
661 while (*ch2 == ' ' || *ch2 == '\t') ch2++;
662
663 while (*ch1 == *ch2)
664 {
665 if (*ch1 == ' ' || *ch1 == '\t' || *ch1 == '\0' || *ch1 == '\n' || *ch1 == '\r') return (1);
666 ch1++;
667 ch2++;
668 }
669
670 if ((*ch1 == ' ' || *ch1 == '\t' || *ch1 == '\0' || *ch1 == '\n' || *ch1 == '\r')
671 && (*ch2 == ' ' || *ch2 == '\t' || *ch2 == '\0' || *ch2 == '\n' || *ch2 == '\r'))
672 return (1);
673 return 0;
674 }
675
676 /* Parses a number in a character string.
677 This routine will detect numbers of the form:
678 nnnn
679 -nnnn
680
681 Args: ch - pointer to the number, in character form.
682 val - integer value returned
683 return value - address of 1st character past the
684 number parsed. NULL if no number found
685 at pointer ch or if the number is an
686 invalid format.
687 */
688
689 char *
intprs(char * ch,int * val)690 intprs(char *ch, int *val)
691 {
692
693 int nflag, flag;
694
695 nflag = 0;
696 if (*ch == '-')
697 {
698 nflag = 1;
699 ch++;
700 }
701 else if (*ch == '+')
702 ch++;
703
704 *val = 0;
705 flag = 1;
706
707 while (*ch >= '0' && *ch <= '9')
708 {
709 *val = *val * 10 + (int) (*ch - '0');
710 flag = 0;
711 ch++;
712 }
713
714 if (flag) return (NULL);
715
716 if (nflag) *val = -1 * *val;
717 return (ch);
718 }
719
720 char *
longprs(char * ch,long * val)721 longprs(char *ch, long *val)
722 {
723
724 int nflag, flag;
725
726 nflag = 0;
727 if (*ch == '-')
728 {
729 nflag = 1;
730 ch++;
731 }
732 else if (*ch == '+')
733 ch++;
734
735 *val = 0;
736 flag = 1;
737
738 while (*ch >= '0' && *ch <= '9')
739 {
740 *val = *val * 10 + (int) (*ch - '0');
741 flag = 0;
742 ch++;
743 }
744
745 if (flag) return (NULL);
746
747 if (nflag) *val = -1 * *val;
748 return (ch);
749 }
750
751 /* Moves a pointer to the start of the next blank-delimited word
752 in a string. If not found, NULL is returned. */
753
754 char *
nxtwrd(char * ch)755 nxtwrd(char *ch)
756 {
757
758 while (*ch != ' ' && *ch != '\t')
759 { /* Skip 1st word */
760 if (*ch == '\0' || *ch == '\n' || *ch == '\r') return (NULL);
761 ch++;
762 }
763 while (*ch == ' ' || *ch == '\t') ch++; /* Find next word */
764 if (*ch == '\0' || *ch == '\n' || *ch == '\r') return (NULL);
765 return (ch);
766 }
767
768 /* Copies a string of a specified length, or when \0 or \n is hit.
769 Trailing blanks are removed, and the output string is terminated
770 with '\0'. */
771
772 void
getstr(char * ch1,char * ch2,int len)773 getstr(char *ch1, char *ch2, int len)
774 {
775 char *ch;
776
777 ch = ch1;
778 while (len > 0 && *ch2 != '\n' && *ch2 != '\0')
779 {
780 *ch1 = *ch2;
781 len--;
782 ch1++;
783 ch2++;
784 }
785 ch1--;
786 while (ch1 >= ch && *ch1 == ' ') ch1--;
787 ch1++;
788 *ch1 = '\0';
789 }
790
791 /* Copies a word of a specified length, or when \0 or \n or \r or ' ' is
792 encountered. The word is terminated with '\0'. ch2 is src, ch1 is dest */
793
794 void
getwrd(char * ch1,char * ch2,int len)795 getwrd(char *ch1, char *ch2, int len)
796 {
797
798 while (len > 0 && *ch2 != '\n' && *ch2 != '\0' && *ch2 != '\r' && *ch2 != ' ' && *ch2 != '\t')
799 {
800 *ch1 = *ch2;
801 len--;
802 ch1++;
803 ch2++;
804 }
805 *ch1 = '\0';
806 }
807
808 /* Determines word length up to next delimiter */
809
810 gaint
wrdlen(char * ch2)811 wrdlen(char *ch2)
812 {
813 gaint len = 0;
814 while (*ch2 != '\n' && *ch2 != '\0' && *ch2 != ' ' && *ch2 != '\t')
815 {
816 len++;
817 ch2++;
818 }
819 return len;
820 }
821
822 /* Converts strings to double */
823 char *
getdbl(char * ch,double * val)824 getdbl(char *ch, double *val)
825 {
826 char *pos = NULL;
827 double res = ch ? strtod(ch, &pos) : 0.0;
828 if (pos == ch)
829 {
830 return NULL;
831 }
832 else
833 {
834 *val = res;
835 return pos;
836 }
837 }
838
839 /* Converts strings to double */
840 char *
getflt(char * ch,float * val)841 getflt(char *ch, float *val)
842 {
843 char *pos = NULL;
844 *val = ch ? (float) strtod(ch, &pos) : 0.0f;
845 if (pos == ch)
846 {
847 return NULL;
848 }
849 else
850 {
851 return pos;
852 }
853 }
854
855 /* Expand file names prefixed with '^' from data descriptor
856 files */
857
858 void
fnmexp(char * out,char * in1,char * in2)859 fnmexp(char *out, char *in1, char *in2)
860 {
861 char *pos, *ch, envv[20], *envr, CR = 13;
862 int i, j;
863
864 if (*in1 == '$')
865 {
866 in1++;
867 i = 0;
868 while (*in1 != '/' && *in1 != '\0' && i < 16)
869 {
870 envv[i] = *in1;
871 i++;
872 in1++;
873 }
874 envv[i] = '\0';
875 envr = getenv(envv);
876 if (envr)
877 {
878 i = 0;
879 j = 0;
880 while (*(envr + j))
881 {
882 *(out + i) = *(envr + j);
883 i++;
884 j++;
885 }
886 /* handle CR for descriptor files created under MS Windows */
887 while (*in1 != '\0' && *in1 != ' ' && *in1 != '\n' && *in1 != CR)
888 {
889 *(out + i) = *in1;
890 i++;
891 in1++;
892 }
893 *(out + i) = '\0';
894 }
895 return;
896 }
897 ch = in2;
898 pos = NULL;
899 while (*ch != '\0' && *ch != ' ' && *ch != '\n')
900 {
901 if (*ch == '/') pos = ch;
902 ch++;
903 }
904 if (pos) pos++;
905 while (pos != NULL && in2 < pos)
906 {
907 *out = *in2;
908 out++;
909 in2++;
910 }
911 in1++;
912 while (*in1 != '\0' && *in1 != ' ' && *in1 != '\n' && *in1 != CR)
913 {
914 *out = *in1;
915 out++;
916 in1++;
917 }
918 *out = '\0';
919 }
920
921 /* Linear conversion routine for dimension conversions. */
922
923 gadouble
liconv(gadouble * vals,gadouble v)924 liconv(gadouble *vals, gadouble v)
925 {
926 return ((*vals * v) + *(vals + 1));
927 }
928
929 /* Non-linear scaling routine for discrete levels. Linear interp
930 between levels. Scaling beyond upper and lower bounds is
931 linear based on the last and first grid spacing, respectively.
932 In each case a pointer to a list of values is provided. The
933 list contains in its first element the number of values
934 in the list. */
935
936 /* Convert a grid value to a world coordinate value.
937 This operation needs to be efficient, since it gets done
938 very often. */
939
940 gadouble
gr2lev(gadouble * vals,gadouble gr)941 gr2lev(gadouble *vals, gadouble gr)
942 {
943 gaint i;
944 if (gr < 1.0) return (*(vals + 1) + (1.0 - gr) * (*(vals + 1) - *(vals + 2)));
945 if (gr > *vals)
946 {
947 i = (gaint)(*vals + 0.1);
948 return (*(vals + i) + (gr - *vals) * (*(vals + i) - *(vals + i - 1)));
949 }
950 i = (gaint) gr;
951 return (*(vals + i) + ((gr - (gadouble) i) * (*(vals + i + 1) - *(vals + i))));
952 }
953
954 /* Convert from world coordinate value to grid value. This operation
955 is not set up to be efficient, under the assumption that it won't
956 get done all that often. */
957
958 gadouble
lev2gr(gadouble * vals,gadouble lev)959 lev2gr(gadouble *vals, gadouble lev)
960 {
961 gaint i, num;
962 gadouble gr;
963 num = (gaint)(*vals + 0.1);
964 for (i = 1; i < num; i++)
965 {
966 if ((lev >= *(vals + i) && lev <= *(vals + i + 1)) || (lev <= *(vals + i) && lev >= *(vals + i + 1)))
967 {
968 gr = (gadouble) i + (lev - *(vals + i)) / (*(vals + i + 1) - *(vals + i));
969 return (gr);
970 }
971 }
972 if (*(vals + 1) < *(vals + num))
973 {
974 if (lev < *(vals + 1))
975 {
976 gr = 1.0 + ((lev - *(vals + 1)) / (*(vals + 2) - *(vals + 1)));
977 return (gr);
978 }
979 gr = (gadouble) i + ((lev - *(vals + i)) / (*(vals + i) - *(vals + i - 1)));
980 return (gr);
981 }
982 else
983 {
984 if (lev > *(vals + 1))
985 {
986 gr = 1.0 + ((lev - *(vals + 1)) / (*(vals + 2) - *(vals + 1)));
987 return (gr);
988 }
989 gr = (gadouble) i + ((lev - *(vals + i)) / (*(vals + i) - *(vals + i - 1)));
990 return (gr);
991 }
992 }
993
994 /* Process linear scaling args */
995
996 gaint
deflin(char * ch,dsets_t * pfi,gaint dim,gaint flag)997 deflin(char *ch, dsets_t *pfi, gaint dim, gaint flag)
998 {
999 gadouble *vals, v1, v2;
1000
1001 vals = (gadouble *) galloc(sizeof(gadouble) * 6, "vals1");
1002 if (vals == NULL) return -1;
1003
1004 if ((ch = nxtwrd(ch)) == NULL) goto err1;
1005 if (getdbl(ch, &v1) == NULL) goto err1;
1006 if (flag)
1007 v2 = 1.0;
1008 else
1009 {
1010 if ((ch = nxtwrd(ch)) == NULL) goto err2;
1011 if (getdbl(ch, &v2) == NULL) goto err2;
1012 }
1013 if (dim != 3 && v2 <= 0.0) goto err2;
1014 *(vals) = v2;
1015 *(vals + 1) = v1 - v2;
1016 *(vals + 2) = -999.9;
1017 pfi->grvals[dim] = vals;
1018 *(vals + 4) = -1.0 * ((v1 - v2) / v2);
1019 *(vals + 3) = 1.0 / v2;
1020 *(vals + 5) = -999.9;
1021 pfi->abvals[dim] = vals + 3;
1022 pfi->ab2gr[dim] = liconv;
1023 pfi->gr2ab[dim] = liconv;
1024 pfi->linear[dim] = 1;
1025 return 0;
1026
1027 err1:
1028 gaprnt(0, "Open Error: Missing or invalid dimension");
1029 gaprnt(0, " starting value\n");
1030 gree(vals, "f178");
1031 return (1);
1032
1033 err2:
1034 gaprnt(0, "Open Error: Missing or invalid dimension");
1035 gaprnt(0, " increment value\n");
1036 gree(vals, "179");
1037 return (1);
1038 }
1039
1040 /* Process levels values in def record */
1041 /* Return codes: -1 is memory allocation error, 1 is other error */
1042
1043 gaint
deflev(char * ch,char * rec,dsets_t * pfi,gaint dim)1044 deflev(char *ch, char *rec, dsets_t *pfi, gaint dim)
1045 {
1046 gadouble *vvs, *vals, v1;
1047 gaint i;
1048
1049 if (pfi->dnum[dim] == 1)
1050 {
1051 i = deflin(ch, pfi, dim, 1);
1052 return (i);
1053 }
1054
1055 vals = (gadouble *) galloc((pfi->dnum[dim] + 5) * sizeof(gadouble), "vals2");
1056 if (vals == NULL) return -1;
1057
1058 vvs = vals;
1059 *vvs = (gadouble) pfi->dnum[dim];
1060 vvs++;
1061 for (i = 0; i < pfi->dnum[dim]; i++)
1062 {
1063 if ((ch = nxtwrd(ch)) == NULL)
1064 {
1065 if (fgets(rec, 256, descr) == NULL) goto err2;
1066 ch = rec;
1067 while (*ch == ' ' || *ch == '\t') ch++;
1068 if (*ch == '\0' || *ch == '\n') goto err3;
1069 }
1070 if (getdbl(ch, &v1) == NULL) goto err1;
1071 *vvs = v1;
1072 vvs++;
1073 }
1074 *vvs = -999.9;
1075 pfi->abvals[dim] = vals;
1076 pfi->grvals[dim] = vals;
1077 pfi->ab2gr[dim] = lev2gr;
1078 pfi->gr2ab[dim] = gr2lev;
1079 pfi->linear[dim] = 0;
1080 return 0;
1081
1082 err1:
1083 gaprnt(0, "Open Error: Invalid value in LEVELS data\n");
1084 gree(vals, "f180");
1085 return (1);
1086
1087 err2:
1088 gaprnt(0, "Open Error: Unexpected EOF reading descriptor file\n");
1089 gaprnt(0, " EOF occurred reading LEVELS values\n");
1090 gree(vals, "f181");
1091 return (1);
1092
1093 err3:
1094 gaprnt(0, "Open Error: Blank Record found in LEVELS data\n");
1095 gree(vals, "f182");
1096 return (1);
1097 }
1098
1099 /* handle var name of the form longnm=>abbrv
1100 or just the abbrv with no long name */
1101
1102 gaint
getvnm(struct gavar * pvar,char * mrec)1103 getvnm(struct gavar *pvar, char *mrec)
1104 {
1105 gaint ib, i, j, k, len, flag;
1106
1107 ib = 0;
1108 while (*(mrec + ib) == ' ') ib++;
1109
1110 if (*(mrec + ib) == '\0' || *(mrec + ib) == '\n') return (1);
1111
1112 /* Scan for the '=>' string */
1113 len = 0;
1114 i = ib;
1115 flag = 0;
1116
1117 while (1)
1118 {
1119 if (*(mrec + i) == ' ' || *(mrec + i) == '\0' || *(mrec + i) == '\n') break;
1120 if (*(mrec + i) == '=' && *(mrec + i + 1) == '>')
1121 {
1122 flag = 1;
1123 break;
1124 }
1125 len++;
1126 i++;
1127 }
1128
1129 if (flag)
1130 {
1131 for (j = ib; j < i; j++)
1132 {
1133 k = j - ib;
1134 pvar->longnm[k] = *(mrec + j);
1135 /* substitute ~ for spaces in longname */
1136 if (pvar->longnm[k] == '~') pvar->longnm[k] = ' ';
1137 }
1138 pvar->longnm[len] = '\0';
1139 i += 2;
1140 }
1141 else
1142 {
1143 i = 0;
1144 pvar->longnm[0] = '\0';
1145 }
1146
1147 if (*(mrec + i) == '\n' || *(mrec + i) == '\0') return (1);
1148
1149 getwrd(pvar->abbrv, mrec + i, 15);
1150 lowcas(pvar->abbrv);
1151
1152 /* Check if 1st character is lower-case alphabetic */
1153 if (islower(*(pvar->abbrv)))
1154 return (0);
1155 else
1156 return (1);
1157 }
1158
1159 /* Given a file name template and a dt structure, fill in to get the file name
1160 */
1161
1162 char *
gafndt(char * fn,struct dt * dtim,struct dt * dtimi,gadouble * vals,struct gachsub * pch1st,struct gaens * ens1st,gaint t,gaint e,gaint * flag)1163 gafndt(char *fn, struct dt *dtim, struct dt *dtimi, gadouble *vals, struct gachsub *pch1st, struct gaens *ens1st, gaint t, gaint e,
1164 gaint *flag)
1165 {
1166 struct gachsub *pchsub;
1167 struct gaens *ens;
1168 gaint len, olen, iv, tdif, i, tused, eused;
1169 char *fnout, *in, *out, *work, *in2, *out2;
1170 (void) vals;
1171 tused = eused = 0;
1172 olen = 0;
1173 while (*(fn + olen)) olen++;
1174 olen += 5;
1175 fnout = (char *) galloc(olen, "fnout");
1176 if (fnout == NULL) return (NULL);
1177
1178 in = fn;
1179 out = fnout;
1180
1181 while (*in)
1182 {
1183 pchsub = pch1st;
1184 ens = ens1st;
1185 /* handle template strings for initial time */
1186 if (*in == '%' && *(in + 1) == 'i')
1187 {
1188 tused = 1;
1189 if (*(in + 2) == 'x' && *(in + 3) == '1')
1190 {
1191 sprintf(out, "%i", dtimi->yr / 10);
1192 while (*out) out++;
1193 in += 4;
1194 }
1195 else if (*(in + 2) == 'x' && *(in + 3) == '3')
1196 {
1197 sprintf(out, "%03i", dtimi->yr / 10);
1198 out += 3;
1199 in += 4;
1200 }
1201 else if (*(in + 2) == 'y' && *(in + 3) == '2')
1202 {
1203 iv = dtimi->yr / 100;
1204 iv = dtimi->yr - iv * 100;
1205 sprintf(out, "%02i", iv);
1206 out += 2;
1207 in += 4;
1208 }
1209 else if (*(in + 2) == 'y' && *(in + 3) == '4')
1210 {
1211 sprintf(out, "%04i", dtimi->yr);
1212 out += 4;
1213 in += 4;
1214 }
1215 else if (*(in + 2) == 'm' && *(in + 3) == '1')
1216 {
1217 sprintf(out, "%i", dtimi->mo);
1218 while (*out) out++;
1219 in += 4;
1220 }
1221 else if (*(in + 2) == 'm' && *(in + 3) == '2')
1222 {
1223 sprintf(out, "%02i", dtimi->mo);
1224 out += 2;
1225 in += 4;
1226 }
1227 else if (*(in + 2) == 'm' && *(in + 3) == 'h')
1228 {
1229 if (dtimi->dy < 16)
1230 *out = 'a';
1231 else
1232 *out = 'b';
1233 out += 1;
1234 in += 4;
1235 }
1236 else if (*(in + 2) == 'm' && *(in + 3) == 'H')
1237 {
1238 if (dtimi->dy < 16)
1239 *out = 'A';
1240 else
1241 *out = 'B';
1242 out += 1;
1243 in += 4;
1244 }
1245 else if (*(in + 2) == 'm' && *(in + 3) == 'c')
1246 {
1247 *out = *(mons[dtimi->mo - 1]);
1248 *(out + 1) = *(mons[dtimi->mo - 1] + 1);
1249 *(out + 2) = *(mons[dtimi->mo - 1] + 2);
1250 out += 3;
1251 in += 4;
1252 }
1253 else if (*(in + 2) == 'd' && *(in + 3) == '1')
1254 {
1255 sprintf(out, "%i", dtimi->dy);
1256 while (*out) out++;
1257 in += 4;
1258 }
1259 else if (*(in + 2) == 'd' && *(in + 3) == '2')
1260 {
1261 sprintf(out, "%02i", dtimi->dy);
1262 out += 2;
1263 in += 4;
1264 }
1265 else if (*(in + 2) == 'h' && *(in + 3) == '1')
1266 {
1267 sprintf(out, "%i", dtimi->hr);
1268 while (*out) out++;
1269 in += 4;
1270 }
1271 else if (*(in + 2) == 'h' && *(in + 3) == '2')
1272 {
1273 sprintf(out, "%02i", dtimi->hr);
1274 out += 2;
1275 in += 4;
1276 }
1277 else if (*(in + 2) == 'h' && *(in + 3) == '3')
1278 {
1279 sprintf(out, "%03i", dtimi->hr);
1280 out += 3;
1281 in += 4;
1282 }
1283 else if (*(in + 2) == 'n' && *(in + 3) == '2')
1284 {
1285 sprintf(out, "%02i", dtimi->mn);
1286 out += 2;
1287 in += 4;
1288 }
1289 else
1290 {
1291 *out = *in;
1292 in++;
1293 out++;
1294 }
1295 }
1296 /* handle template strings for any time */
1297 else if (*in == '%' && *(in + 1) == 'x' && *(in + 2) == '1')
1298 { /* x: decades */
1299 tused = 1;
1300 sprintf(out, "%i", dtim->yr / 10);
1301 while (*out) out++;
1302 in += 3;
1303 }
1304 else if (*in == '%' && *(in + 1) == 'x' && *(in + 2) == '3')
1305 {
1306 tused = 1;
1307 sprintf(out, "%03i", dtim->yr / 10);
1308 out += 3;
1309 in += 3;
1310 }
1311 else if (*in == '%' && *(in + 1) == 'y' && *(in + 2) == '2')
1312 {
1313 tused = 1;
1314 iv = dtim->yr / 100;
1315 iv = dtim->yr - iv * 100;
1316 sprintf(out, "%02i", iv);
1317 out += 2;
1318 in += 3;
1319 }
1320 else if (*in == '%' && *(in + 1) == 'y' && *(in + 2) == '4')
1321 {
1322 tused = 1;
1323 sprintf(out, "%04i", dtim->yr);
1324 out += 4;
1325 in += 3;
1326 }
1327 else if (*in == '%' && *(in + 1) == 'm' && *(in + 2) == '1')
1328 {
1329 tused = 1;
1330 sprintf(out, "%i", dtim->mo);
1331 while (*out) out++;
1332 in += 3;
1333 }
1334 else if (*in == '%' && *(in + 1) == 'm' && *(in + 2) == '2')
1335 {
1336 tused = 1;
1337 sprintf(out, "%02i", dtim->mo);
1338 out += 2;
1339 in += 3;
1340 }
1341 else if (*in == '%' && *(in + 1) == 'm' && *(in + 2) == 'h')
1342 {
1343 tused = 1;
1344 if (dtim->dy < 16)
1345 *out = 'a';
1346 else
1347 *out = 'b';
1348 out += 1;
1349 in += 3;
1350 }
1351 else if (*in == '%' && *(in + 1) == 'm' && *(in + 2) == 'H')
1352 {
1353 tused = 1;
1354 if (dtim->dy < 16)
1355 *out = 'A';
1356 else
1357 *out = 'B';
1358 out += 1;
1359 in += 3;
1360 }
1361 else if (*in == '%' && *(in + 1) == 'm' && *(in + 2) == 'c')
1362 {
1363 tused = 1;
1364 *out = *(mons[dtim->mo - 1]);
1365 *(out + 1) = *(mons[dtim->mo - 1] + 1);
1366 *(out + 2) = *(mons[dtim->mo - 1] + 2);
1367 out += 3;
1368 in += 3;
1369 }
1370 else if (*in == '%' && *(in + 1) == 'd' && *(in + 2) == '1')
1371 {
1372 tused = 1;
1373 sprintf(out, "%i", dtim->dy);
1374 while (*out) out++;
1375 in += 3;
1376 }
1377 else if (*in == '%' && *(in + 1) == 'd' && *(in + 2) == '2')
1378 {
1379 tused = 1;
1380 sprintf(out, "%02i", dtim->dy);
1381 out += 2;
1382 in += 3;
1383 }
1384 else if (*in == '%' && *(in + 1) == 'h' && *(in + 2) == '1')
1385 {
1386 tused = 1;
1387 sprintf(out, "%i", dtim->hr);
1388 while (*out) out++;
1389 in += 3;
1390 }
1391 else if (*in == '%' && *(in + 1) == 'h' && *(in + 2) == '2')
1392 {
1393 tused = 1;
1394 sprintf(out, "%02i", dtim->hr);
1395 out += 2;
1396 in += 3;
1397 }
1398 else if (*in == '%' && *(in + 1) == 'h' && *(in + 2) == '3')
1399 {
1400 tused = 1;
1401 sprintf(out, "%03i", dtim->hr);
1402 out += 3;
1403 in += 3;
1404 }
1405 else if (*in == '%' && *(in + 1) == 'n' && *(in + 2) == '2')
1406 {
1407 tused = 1;
1408 sprintf(out, "%02i", dtim->mn);
1409 out += 2;
1410 in += 3;
1411 }
1412 /* forecast times */
1413 else if (*in == '%' && *(in + 1) == 'f' && *(in + 2) == '2')
1414 {
1415 tused = 1;
1416 tdif = timdif(dtimi, dtim);
1417 tdif = (tdif + 30) / 60;
1418 if (tdif < 99)
1419 sprintf(out, "%02i", tdif);
1420 else
1421 sprintf(out, "%i", tdif);
1422 while (*out) out++;
1423 in += 3;
1424 }
1425 else if (*in == '%' && *(in + 1) == 'f' && *(in + 2) == '3')
1426 {
1427 tused = 1;
1428 tdif = timdif(dtimi, dtim);
1429 tdif = (tdif + 30) / 60;
1430 if (tdif < 999)
1431 sprintf(out, "%03i", tdif);
1432 else
1433 sprintf(out, "%i", tdif);
1434 while (*out) out++;
1435 in += 3;
1436 }
1437 /* string substitution */
1438 else if (*in == '%' && *(in + 1) == 'c' && *(in + 2) == 'h')
1439 {
1440 tused = 1;
1441 while (pchsub)
1442 {
1443 if (t >= pchsub->t1 && (pchsub->t2 == -99 || t <= pchsub->t2))
1444 {
1445 len = wrdlen(pchsub->ch); /* Reallocate output string */
1446 olen += len;
1447 work = (char *) galloc(olen, "work");
1448 if (work == NULL)
1449 {
1450 gree(fnout, "f240");
1451 return (NULL);
1452 }
1453 in2 = fnout;
1454 out2 = work;
1455 while (in2 != out)
1456 {
1457 *out2 = *in2;
1458 in2++;
1459 out2++;
1460 }
1461 gree(fnout, "f241");
1462 fnout = work;
1463 out = out2;
1464 getwrd(out, pchsub->ch, len);
1465 out += len;
1466 break;
1467 }
1468 pchsub = pchsub->forw;
1469 }
1470 in += 3;
1471 }
1472 /* ensemble name substitution */
1473 else if (*in == '%' && *(in + 1) == 'e')
1474 {
1475 eused = 1;
1476 if (ens == NULL)
1477 {
1478 gree(fnout, "f242");
1479 return (NULL);
1480 }
1481 else
1482 {
1483 /* advance through array of ensemble structures, till we reach
1484 * ensemble 'e' */
1485 i = 1;
1486 while (i != e)
1487 {
1488 i++;
1489 ens++;
1490 }
1491 len = strlen(ens->name);
1492 if (len < 1)
1493 {
1494 gree(fnout, "f243");
1495 return (NULL);
1496 }
1497 olen += len;
1498 work = (char *) galloc(olen, "work2"); /* Reallocate output string */
1499 if (work == NULL)
1500 {
1501 gree(fnout, "f244");
1502 return (NULL);
1503 }
1504 in2 = fnout; /* copy the string we've got so far */
1505 out2 = work;
1506 while (in2 != out)
1507 {
1508 *out2 = *in2;
1509 in2++;
1510 out2++;
1511 }
1512 gree(fnout, "f245");
1513 fnout = work;
1514 out = out2;
1515 getwrd(out, ens->name, len);
1516 out += len;
1517 }
1518 in += 2;
1519 }
1520 else
1521 {
1522 *out = *in;
1523 in++;
1524 out++;
1525 }
1526 }
1527 *out = '\0';
1528 if (eused == 1 && tused == 1)
1529 {
1530 *flag = 3; /* templating on E and T */
1531 }
1532 else if (eused == 1 && tused == 0)
1533 {
1534 *flag = 2; /* templating only on E */
1535 }
1536 else if (eused == 0 && tused == 1)
1537 {
1538 *flag = 1; /* templating only on T */
1539 }
1540 else
1541 {
1542 *flag = 0; /* no templating */
1543 }
1544 return (fnout);
1545 }
1546
1547 #undef IsBigendian
1548 #define IsBigendian() (u_byteorder.c[sizeof(long) - 1])
1549
1550 int
read_gradsdes(char * filename,dsets_t * pfi)1551 read_gradsdes(char *filename, dsets_t *pfi)
1552 {
1553 /* IsBigendian returns 1 for big endian byte order */
1554 static union
1555 {
1556 unsigned long l;
1557 unsigned char c[sizeof(long)];
1558 } u_byteorder = { 1 };
1559 struct gavar *pvar;
1560 struct gaens *ens;
1561 struct dt tdef, tdefe, tdefi, dt1, dt2;
1562 struct gachsub *pchsub;
1563 int status = 0;
1564 int reclen;
1565 int ichar;
1566 char rec[MAX_RECLEN], mrec[MAX_RECLEN];
1567 char *ch, *pos;
1568 int i, j, ii, jj;
1569 off_t levs, acum, recacm;
1570 // gaint acumstride=0;
1571 gaint hdrb, trlb;
1572 gaint size = 0, rc, len, flag, tim1, tim2;
1573 gaint e, t;
1574 int BYTEORDER = IsBigendian();
1575 gadouble *vals;
1576 gadouble v1, v2, temp;
1577 int err = 0;
1578
1579 hdrb = 0;
1580 trlb = 0;
1581
1582 /* Try to open descriptor file */
1583 descr = fopen(filename, "r");
1584 if (descr == NULL)
1585 {
1586 fprintf(stderr, "fopen failed on %s", filename);
1587 return -1;
1588 }
1589
1590 /* Copy descriptor file name into gafile structure */
1591 getwrd(pfi->dnam, filename, MAX_NAMELEN);
1592
1593 /* Parse the data descriptor file */
1594 while (fgets(rec, MAX_RECLEN, descr) != NULL)
1595 {
1596
1597 /* Remove any leading blanks from rec */
1598 reclen = (int) strlen(rec);
1599 jj = 0;
1600 while (jj < reclen && rec[0] == ' ')
1601 {
1602 for (ii = 0; ii < reclen; ii++) rec[ii] = rec[ii + 1];
1603 jj++;
1604 }
1605 /* replace newline with null at end of record */
1606 for (ichar = (int) strlen(rec) - 1; ichar >= 0; --ichar)
1607 {
1608 if (rec[ichar] == '\n')
1609 {
1610 rec[ichar] = '\0';
1611 break;
1612 }
1613 }
1614 /* Keep mixed case and lower case versions of rec handy */
1615 strcpy(mrec, rec);
1616 lowcas(rec);
1617
1618 if (!isalnum(mrec[0]))
1619 {
1620 /* check if comment contains attribute metadata */
1621 /*
1622 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0))
1623 {
1624 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
1625 }
1626 */
1627 }
1628 else if (cmpwrd("byteswapped", rec))
1629 {
1630 pfi->bswap = 1;
1631 }
1632 else if (cmpwrd("fileheader", rec))
1633 {
1634 if ((ch = nxtwrd(rec)) == NULL)
1635 {
1636 gaprnt(1, "Description file warning: Missing fileheader length\n");
1637 }
1638 else
1639 {
1640 ch = longprs(ch, &(pfi->fhdr));
1641 if (ch == NULL)
1642 {
1643 gaprnt(1, "Fileheader record invalid\n");
1644 pfi->fhdr = 0;
1645 }
1646 }
1647 }
1648 else if (cmpwrd("xyheader", rec))
1649 {
1650 if ((ch = nxtwrd(rec)) == NULL)
1651 {
1652 gaprnt(1, "Description file warning: Missing xy grid header length\n");
1653 }
1654 else
1655 {
1656 ch = longprs(ch, &(pfi->xyhdr));
1657 if (ch == NULL)
1658 {
1659 gaprnt(1, "xy grid header length invalid\n");
1660 pfi->xyhdr = 0;
1661 }
1662 else
1663 {
1664 pfi->xyhdr = pfi->xyhdr / 4;
1665 }
1666 }
1667 }
1668 else if (cmpwrd("format", rec) || cmpwrd("options", rec))
1669 {
1670 if ((ch = nxtwrd(rec)) == NULL)
1671 {
1672 gaprnt(1, "Description file warning: Missing options keyword\n");
1673 }
1674 else
1675 {
1676 while (ch != NULL)
1677 {
1678 if (cmpwrd("sequential", ch))
1679 pfi->seqflg = 1;
1680 else if (cmpwrd("yrev", ch))
1681 pfi->yrflg = 1;
1682 else if (cmpwrd("zrev", ch))
1683 pfi->zrflg = 1;
1684 else if (cmpwrd("template", ch))
1685 pfi->tmplat = 1;
1686 else if (cmpwrd("flt64", ch))
1687 pfi->flt64 = 1;
1688 else if (cmpwrd("byteswapped", ch))
1689 pfi->bswap = 1;
1690 #if GRIB2
1691 else if (cmpwrd("pascals", ch))
1692 pfi->pa2mb = 1;
1693 #endif
1694 else if (cmpwrd("365_day_calendar", ch))
1695 {
1696 pfi->calendar = 1;
1697 }
1698 else if (cmpwrd("big_endian", ch))
1699 {
1700 if (!BYTEORDER) pfi->bswap = 1;
1701 }
1702 else if (cmpwrd("little_endian", ch))
1703 {
1704 if (BYTEORDER) pfi->bswap = 1;
1705 }
1706 else
1707 {
1708 gaprnt(0, "Open Error: Data file type invalid\n");
1709 goto err9;
1710 }
1711 ch = nxtwrd(ch);
1712 }
1713 }
1714 }
1715 else if (cmpwrd("trailerbytes", rec))
1716 {
1717 if ((ch = nxtwrd(rec)) == NULL)
1718 {
1719 gaprnt(1, "Trailerbytes record invalid\n");
1720 }
1721 else
1722 {
1723 ch = intprs(ch, &trlb);
1724 if (ch == NULL)
1725 {
1726 gaprnt(1, "Trailerbytes record invalid\n");
1727 trlb = 0;
1728 }
1729 else
1730 {
1731 trlb = trlb / 4;
1732 }
1733 }
1734 }
1735 else if (cmpwrd("headerbytes", rec) || cmpwrd("theader", rec))
1736 {
1737 if ((ch = nxtwrd(rec)) == NULL)
1738 {
1739 gaprnt(1, "headerbytes/theader record invalid\n");
1740 }
1741 else
1742 {
1743 ch = intprs(ch, &hdrb);
1744 if (ch == NULL)
1745 {
1746 gaprnt(1, "headerbytes/theader record invalid\n");
1747 hdrb = 0;
1748 }
1749 else
1750 {
1751 hdrb = hdrb / 4;
1752 }
1753 }
1754 }
1755 /* Handle the chsub records. time1, time2, then a string, multiple times
1756 */
1757 else if (cmpwrd("chsub", rec))
1758 {
1759 /* point to first block in chain */
1760 pchsub = pfi->pchsub1;
1761 if (pchsub != NULL)
1762 {
1763 while (pchsub->forw != NULL)
1764 {
1765 pchsub = pchsub->forw; /* advance to end of chain */
1766 }
1767 }
1768 flag = 0;
1769 ch = mrec;
1770 while (1)
1771 {
1772 if ((ch = nxtwrd(ch)) == NULL) break;
1773 flag = 1;
1774 if ((ch = intprs(ch, &tim1)) == NULL) break;
1775 if ((ch = nxtwrd(ch)) == NULL) break;
1776 if (*ch == '*' && (*(ch + 1) == ' ' || *(ch + 1) == '\t'))
1777 tim2 = -99;
1778 else if ((ch = intprs(ch, &tim2)) == NULL)
1779 break;
1780 if ((ch = nxtwrd(ch)) == NULL) break;
1781 flag = 0;
1782 if (pchsub)
1783 { /* chain exists */
1784 pchsub->forw = (struct gachsub *) galloc(sizeof(struct gachsub), "chsubnew");
1785 if (pchsub->forw == NULL)
1786 {
1787 gaprnt(0, "Open Error: memory allocation failed for pchsub\n");
1788 goto err8;
1789 }
1790 pchsub = pchsub->forw;
1791 pchsub->forw = NULL;
1792 }
1793 else
1794 { /* start a new chain */
1795 pfi->pchsub1 = (struct gachsub *) galloc(sizeof(struct gachsub), "chsub1");
1796 if (pfi->pchsub1 == NULL)
1797 {
1798 gaprnt(0, "Open Error: memory allocation failed for pchsub1\n");
1799 goto err8;
1800 }
1801 pchsub = pfi->pchsub1;
1802 pchsub->forw = NULL;
1803 }
1804 len = wrdlen(ch);
1805 if ((pchsub->ch = (char *) galloc(len + 1, "chsubstr")) == NULL) goto err8;
1806 getwrd(pchsub->ch, ch, len);
1807 pchsub->t1 = tim1;
1808 pchsub->t2 = tim2;
1809 }
1810 if (flag)
1811 {
1812 gaprnt(1, "Description file warning: Invalid chsub record; Ignored\n");
1813 }
1814 }
1815 else if (cmpwrd("title", rec))
1816 {
1817 if ((ch = nxtwrd(mrec)) == NULL)
1818 {
1819 gaprnt(1, "Description file warning: Missing title string\n");
1820 }
1821 else
1822 {
1823 getstr(pfi->title, ch, MAX_NAMELEN);
1824 }
1825 }
1826 else if (cmpwrd("dset", rec))
1827 {
1828 ch = nxtwrd(mrec);
1829 if (ch == NULL)
1830 {
1831 gaprnt(0, "Descriptor File Error: Data file name is missing\n");
1832 goto err9;
1833 }
1834 if (*ch == '^' || *ch == '$')
1835 {
1836 fnmexp(pfi->name, ch, filename);
1837 }
1838 else
1839 {
1840 getwrd(pfi->name, ch, MAX_NAMELEN);
1841 }
1842 }
1843 else if (cmpwrd("undef", rec))
1844 {
1845 ch = nxtwrd(mrec);
1846 if (ch == NULL)
1847 {
1848 gaprnt(0, "Open Error: Missing undef value\n");
1849 goto err9;
1850 }
1851
1852 pos = getdbl(ch, &(pfi->undef));
1853 if (pos == NULL)
1854 {
1855 gaprnt(0, "Open Error: Invalid undef value\n");
1856 goto err9;
1857 }
1858
1859 pfi->ulow = fabs(pfi->undef / EPSILON);
1860 pfi->uhi = pfi->undef + pfi->ulow;
1861 pfi->ulow = pfi->undef - pfi->ulow;
1862 }
1863 else if (cmpwrd("xdef", rec))
1864 {
1865 if (pfi->type == 2) continue;
1866 if ((ch = nxtwrd(rec)) == NULL) goto err1;
1867 if ((pos = intprs(ch, &(pfi->dnum[0]))) == NULL) goto err1;
1868 if (pfi->dnum[0] < 1)
1869 {
1870 sprintf(pout, "Warning: Invalid XDEF syntax in %s -- Changing size of X axis from %d to 1 \n", pfi->dnam,
1871 pfi->dnum[0]);
1872 gaprnt(1, pout);
1873 pfi->dnum[0] = 1;
1874 }
1875 if (*pos != ' ') goto err1;
1876 if ((ch = nxtwrd(ch)) == NULL) goto err2;
1877 if (cmpwrd("linear", ch))
1878 {
1879 rc = deflin(ch, pfi, 0, 0);
1880 if (rc == -1) goto err8;
1881 if (rc) goto err9;
1882 v2 = *(pfi->grvals[0]);
1883 v1 = *(pfi->grvals[0] + 1) + v2;
1884 temp = v1 + ((gadouble)(pfi->dnum[0])) * v2;
1885 temp = temp - 360.0;
1886 if (fabs(temp - v1) < 0.01) pfi->wrap = 1;
1887 }
1888 else if (cmpwrd("levels", ch))
1889 {
1890 rc = deflev(ch, rec, pfi, 0);
1891 if (rc == -1) goto err8;
1892 if (rc) goto err9;
1893 }
1894 else
1895 goto err2;
1896 }
1897 else if (cmpwrd("ydef", rec))
1898 {
1899 if (pfi->type == 2) continue;
1900 if ((ch = nxtwrd(rec)) == NULL) goto err1;
1901 if ((pos = intprs(ch, &(pfi->dnum[1]))) == NULL) goto err1;
1902 if (pfi->dnum[1] < 1)
1903 {
1904 sprintf(pout, "Warning: Invalid YDEF syntax in %s -- Changing size of Y axis from %d to 1 \n", pfi->dnam,
1905 pfi->dnum[1]);
1906 gaprnt(1, pout);
1907 pfi->dnum[1] = 1;
1908 }
1909 if (*pos != ' ') goto err1;
1910 if ((ch = nxtwrd(ch)) == NULL) goto err2;
1911 if (cmpwrd("linear", ch))
1912 {
1913 rc = deflin(ch, pfi, 1, 0);
1914 if (rc == -1) goto err8;
1915 if (rc) goto err9;
1916 }
1917 else if (cmpwrd("levels", ch))
1918 {
1919 rc = deflev(ch, rec, pfi, 1);
1920 if (rc == -1) goto err8;
1921 if (rc) goto err9;
1922 }
1923 }
1924 else if (cmpwrd("zdef", rec))
1925 {
1926 if (pfi->type == 2) continue;
1927 if ((ch = nxtwrd(rec)) == NULL) goto err1;
1928 if ((pos = intprs(ch, &(pfi->dnum[2]))) == NULL) goto err1;
1929 if (pfi->dnum[2] < 1)
1930 {
1931 sprintf(pout, "Warning: Invalid ZDEF syntax in %s -- Changing size of Z axis from %d to 1 \n", pfi->dnam,
1932 pfi->dnum[2]);
1933 gaprnt(1, pout);
1934 pfi->dnum[2] = 1;
1935 }
1936 if (*pos != ' ') goto err1;
1937 if ((ch = nxtwrd(ch)) == NULL) goto err2;
1938 if (cmpwrd("linear", ch))
1939 {
1940 rc = deflin(ch, pfi, 2, 0);
1941 if (rc == -1) goto err8;
1942 if (rc) goto err9;
1943 }
1944 else if (cmpwrd("levels", ch))
1945 {
1946 rc = deflev(ch, rec, pfi, 2);
1947 if (rc == -1) goto err8;
1948 if (rc) goto err9;
1949 }
1950 else
1951 goto err2;
1952 }
1953 else if (cmpwrd("tdef", rec))
1954 {
1955 if ((ch = nxtwrd(rec)) == NULL) goto err1;
1956 if ((pos = intprs(ch, &(pfi->dnum[3]))) == NULL) goto err1;
1957 if (pfi->dnum[3] < 1)
1958 {
1959 sprintf(pout, "Warning: Invalid TDEF syntax in %s -- Changing size of T axis from %d to 1 \n", pfi->dnam,
1960 pfi->dnum[3]);
1961 gaprnt(1, pout);
1962 pfi->dnum[3] = 1;
1963 }
1964 if (*pos != ' ') goto err1;
1965 if ((ch = nxtwrd(ch)) == NULL) goto err2;
1966 if (cmpwrd("linear", ch))
1967 {
1968 if ((ch = nxtwrd(ch)) == NULL) goto err3a_tdef;
1969 tdef.yr = -1000;
1970 tdef.mo = -1000;
1971 tdef.dy = -1000;
1972 if ((pos = adtprs(ch, &tdef, &dt1)) == NULL) goto err3b_tdef;
1973 if (*pos != ' ' || dt1.yr == -1000 || dt1.mo == -1000 || dt1.dy == -1000) goto err3c_tdef;
1974 if ((ch = nxtwrd(ch)) == NULL) goto err4a_tdef;
1975 if ((pos = rdtprs(ch, &dt2)) == NULL) goto err4b_tdef;
1976 v1 = (dt2.yr * 12) + dt2.mo;
1977 v2 = (dt2.dy * 1440) + (dt2.hr * 60) + dt2.mn;
1978 /* check if 0 dt */
1979 if (((int) v1 == 0) && ((int) v2 == 0)) goto err4c_tdef;
1980 if ((vals = (gadouble *) galloc(sizeof(gadouble) * 8, "tvals5")) == NULL) goto err8;
1981 *(vals) = dt1.yr;
1982 *(vals + 1) = dt1.mo;
1983 *(vals + 2) = dt1.dy;
1984 *(vals + 3) = dt1.hr;
1985 *(vals + 4) = dt1.mn;
1986 *(vals + 5) = v1;
1987 *(vals + 6) = v2;
1988 *(vals + 7) = -999.9;
1989 pfi->grvals[3] = vals;
1990 pfi->abvals[3] = vals;
1991 pfi->linear[3] = 1;
1992 }
1993 else
1994 goto err2;
1995 }
1996 else if (cmpwrd("vars", rec))
1997 {
1998 if ((ch = nxtwrd(rec)) == NULL) goto err5;
1999 if ((pos = intprs(ch, &(pfi->vnum))) == NULL) goto err5;
2000 size = pfi->vnum * (sizeof(struct gavar) + 7);
2001 if ((pvar = (struct gavar *) galloc(size, "pvar2")) == NULL) goto err8;
2002 pfi->pvar1 = pvar;
2003 i = 0;
2004 while (i < pfi->vnum)
2005 {
2006 /* initialize variables in the pvar structure */
2007 pvar->offset = 0;
2008 pvar->recoff = 0;
2009 pvar->ncvid = -999;
2010 pvar->sdvid = -999;
2011 pvar->levels = 0;
2012 pvar->dfrm = 0;
2013 pvar->var_t = 0;
2014 pvar->scale = 1;
2015 pvar->add = 0;
2016 pvar->undef = -9.99E33;
2017 pvar->vecpair = -999;
2018 pvar->isu = 0;
2019 pvar->isdvar = 0;
2020 pvar->nvardims = 0;
2021
2022 /* get the complete variable declaration */
2023 if (fgets(rec, 512, descr) == NULL)
2024 {
2025 gaprnt(0, "Open Error: Unexpected EOF reading variables\n");
2026 sprintf(pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
2027 gaprnt(2, pout);
2028 goto retrn;
2029 }
2030 /* remove any leading blanks from rec */
2031 reclen = strlen(rec);
2032 jj = 0;
2033 while (jj < reclen && rec[0] == ' ')
2034 {
2035 for (ii = 0; ii < reclen; ii++) rec[ii] = rec[ii + 1];
2036 jj++;
2037 }
2038 /* replace newline with null at end of record */
2039 for (ichar = strlen(rec) - 1; ichar >= 0; --ichar)
2040 {
2041 if (rec[ichar] == '\n')
2042 {
2043 rec[ichar] = '\0';
2044 break;
2045 }
2046 }
2047 /* Keep mixed case and lower case versions of rec handy */
2048 strcpy(mrec, rec);
2049 lowcas(rec);
2050 /* Allow comments between VARS and ENDVARS */
2051 if (!isalnum(*(mrec)))
2052 {
2053 /* Parse comment if it contains attribute metadata */
2054 /*
2055 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0)) {
2056 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
2057 else continue;
2058 }
2059 else */
2060 continue;
2061 }
2062 if (cmpwrd("endvars", rec))
2063 {
2064 gaprnt(0, "Open Error: Unexpected ENDVARS record\n");
2065 sprintf(pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
2066 gaprnt(2, pout);
2067 goto err9;
2068 }
2069
2070 /* get abbrv and full variable name if there */
2071 if ((getvnm(pvar, mrec)) != 0) goto err6;
2072
2073 /* parse the levels fields */
2074 if ((ch = nxtwrd(rec)) == NULL) goto err6;
2075 /* begin with 8th element of units aray for levels values */
2076 for (j = 0; j < 16; j++) pvar->units[j] = -999;
2077 j = 8;
2078 while (1)
2079 {
2080 if (j == 8)
2081 {
2082 /* first element is num levels */
2083 if ((ch = intprs(ch, &(pvar->levels))) == NULL) goto err6;
2084 }
2085 else
2086 {
2087 /* remaining elements are grib2 level codes */
2088 if ((ch = getdbl(ch, &(pvar->units[j - 1]))) == NULL) goto err6;
2089 }
2090 /* advance through comma-delimited list of levels args */
2091 while (*ch == ' ') ch++;
2092 if (*ch == '\0' || *ch == '\n') goto err6;
2093 if (*ch != ',') break;
2094 ch++;
2095 while (*ch == ',')
2096 {
2097 ch++;
2098 j++;
2099 } /* advance past back to back commas */
2100 while (*ch == ' ') ch++;
2101 if (*ch == '\0' || *ch == '\n') goto err6;
2102 j++;
2103 if (j > 15) goto err6;
2104 }
2105
2106 /* parse the units fields; begin with 0th element for variable
2107 * units */
2108 j = 0;
2109 pvar->nvardims = 0;
2110 while (1)
2111 {
2112 if (*ch == 'x' || *ch == 'y' || *ch == 'z' || *ch == 't' || *ch == 'e')
2113 {
2114 if (*(ch + 1) != ',' && *(ch + 1) != ' ') goto err6;
2115 if (*ch == 'x')
2116 {
2117 pvar->units[j] = -100;
2118 pvar->nvardims++;
2119 }
2120 if (*ch == 'y')
2121 {
2122 pvar->units[j] = -101;
2123 pvar->nvardims++;
2124 }
2125 if (*ch == 'z')
2126 {
2127 pvar->units[j] = -102;
2128 pvar->nvardims++;
2129 }
2130 if (*ch == 't')
2131 {
2132 pvar->units[j] = -103;
2133 pvar->nvardims++;
2134 }
2135 if (*ch == 'e')
2136 {
2137 pvar->units[j] = -104;
2138 pvar->nvardims++;
2139 }
2140 ch++;
2141 }
2142 else
2143 {
2144 if ((ch = getdbl(ch, &(pvar->units[j]))) == NULL) goto err6;
2145 /* no negative array indices for ncflag files */
2146 if ((pfi->ncflg) && (pvar->units[j] < 0)) goto err6;
2147 }
2148 while (*ch == ' ') ch++;
2149 if (*ch == '\0' || *ch == '\n') goto err6;
2150 if (*ch != ',') break;
2151 ch++;
2152 while (*ch == ' ') ch++;
2153 if (*ch == '\0' || *ch == '\n') goto err6;
2154 j++;
2155 if (j > 8) goto err6;
2156 }
2157
2158 /* parse the variable description */
2159 getstr(pvar->varnm, mrec + (ch - rec), 127);
2160
2161 /* var_t is for data files with dimension sequence: X, Y, Z, T, V
2162 */
2163 if (((int) pvar->units[0] == -1) && ((int) pvar->units[1] == 20)) pvar->var_t = 1;
2164
2165 /* non-float data types */
2166 if (((int) pvar->units[0] == -1) && ((int) pvar->units[1] == 40))
2167 {
2168
2169 if ((int) pvar->units[2] == 1) pvar->dfrm = 1;
2170 if ((int) pvar->units[2] == 2)
2171 {
2172 pvar->dfrm = 2;
2173 if ((int) pvar->units[3] == -1) pvar->dfrm = -2;
2174 }
2175 if ((int) pvar->units[2] == 4) pvar->dfrm = 4;
2176 }
2177
2178 i++;
2179 pvar++;
2180 }
2181
2182 /* Get ENDVARS statement and any additional comments */
2183 if (fgets(rec, 512, descr) == NULL)
2184 {
2185 gaprnt(0, "Open Error: Missing ENDVARS statement.\n");
2186 goto retrn;
2187 }
2188 /* Remove any leading blanks from rec */
2189 reclen = strlen(rec);
2190 jj = 0;
2191 while (jj < reclen && rec[0] == ' ')
2192 {
2193 for (ii = 0; ii < reclen; ii++) rec[ii] = rec[ii + 1];
2194 jj++;
2195 }
2196 /* replace newline with null at end of record */
2197 for (ichar = strlen(rec) - 1; ichar >= 0; --ichar)
2198 {
2199 if (rec[ichar] == '\n')
2200 {
2201 rec[ichar] = '\0';
2202 break;
2203 }
2204 }
2205 /* Keep mixed case and lower case versions handy */
2206 strcpy(mrec, rec);
2207 lowcas(rec);
2208 while (!cmpwrd("endvars", rec))
2209 {
2210 /* see if it's an attribute comment */
2211 if (!isalnum(*(mrec)))
2212 {
2213 /*
2214 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0))
2215 { if ((ddfattr(mrec,pfi)) == -1) goto retrn;
2216 }
2217 */
2218 }
2219 else
2220 {
2221 sprintf(pout, "Open Error: Looking for \"endvars\", found \"%s\" instead.\n", rec);
2222 gaprnt(0, pout);
2223 goto err9;
2224 }
2225 /* get a new record */
2226 if (fgets(rec, 512, descr) == NULL)
2227 {
2228 gaprnt(0, "Open Error: Missing ENDVARS statement.\n");
2229 goto retrn;
2230 }
2231 /* Remove any leading blanks from new record */
2232 reclen = strlen(rec);
2233 jj = 0;
2234 while (jj < reclen && rec[0] == ' ')
2235 {
2236 for (ii = 0; ii < reclen; ii++) rec[ii] = rec[ii + 1];
2237 jj++;
2238 }
2239 /* replace newline with null at end of record */
2240 for (ichar = strlen(rec) - 1; ichar >= 0; --ichar)
2241 {
2242 if (rec[ichar] == '\n')
2243 {
2244 rec[ichar] = '\0';
2245 break;
2246 }
2247 }
2248 /* Keep mixed case and lower case versions handy */
2249 strcpy(mrec, rec);
2250 lowcas(rec);
2251 }
2252 /* vars block parsed without error */
2253 }
2254 else
2255 {
2256 /* parse error of .ctl file */
2257 gaprnt(0, "Open Error: Unknown keyword in description file\n");
2258 goto err9;
2259 }
2260 }
2261 /* Done scanning!
2262 Check if scanned stuff makes sense, and then set things up correctly */
2263
2264 pfi->ulow = fabs(pfi->undef / EPSILON);
2265 pfi->uhi = pfi->undef + pfi->ulow;
2266 pfi->ulow = pfi->undef - pfi->ulow;
2267
2268 /* If no EDEF entry was found, set up the default values */
2269 if (pfi->ens1 == NULL)
2270 {
2271 pfi->dnum[4] = 1;
2272 /* set up linear scaling */
2273 if ((vals = (gadouble *) galloc(sizeof(gadouble) * 6, "evals3")) == NULL) goto err8;
2274 v1 = v2 = 1;
2275 *(vals + 1) = v1 - v2;
2276 *(vals) = v2;
2277 *(vals + 2) = -999.9;
2278 pfi->grvals[4] = vals;
2279 *(vals + 4) = -1.0 * ((v1 - v2) / v2);
2280 *(vals + 3) = 1.0 / v2;
2281 *(vals + 5) = -999.9;
2282 pfi->abvals[4] = vals + 3;
2283 pfi->ab2gr[4] = liconv;
2284 pfi->gr2ab[4] = liconv;
2285 pfi->linear[4] = 1;
2286 /* Allocate memory and initialize one ensemble structure */
2287 ens = (struct gaens *) galloc(sizeof(struct gaens), "ens5");
2288 if (ens == NULL)
2289 {
2290 gaprnt(0, "Open Error: memory allocation failed for default ens\n");
2291 goto err8;
2292 }
2293 pfi->ens1 = ens;
2294 sprintf(ens->name, "1");
2295 ens->length = pfi->dnum[3];
2296 ens->gt = 1;
2297 gr2t(pfi->grvals[3], 1, &ens->tinit);
2298 for (j = 0; j < 4; j++) ens->grbcode[j] = -999;
2299 }
2300
2301 /* Make sure there are no conflicting options and data types */
2302 pvar = pfi->pvar1;
2303 for (j = 1; j <= pfi->vnum; j++)
2304 {
2305 if ((int) pvar->units[0] == -1 && (int) pvar->units[1] == 20)
2306 {
2307 if (pfi->tmplat)
2308 {
2309 gaprnt(0, "Open Error: Variables with transposed VAR-T dimensions cannot be templated together\n");
2310 err = 1;
2311 }
2312 if (hdrb > 0)
2313 {
2314 gaprnt(0, "Open Error: Variables with transposed VAR-T dimensions are incompatible with time headers\n");
2315 err = 1;
2316 }
2317 if (trlb > 0)
2318 {
2319 gaprnt(0, "Open Error: Variables with transposed VAR-T dimensions are incompatible with TRAILERBYTES\n");
2320 err = 1;
2321 }
2322 }
2323 pvar++;
2324 }
2325 if (err) goto retrn;
2326
2327 /* Figure out locations of variables within a time group */
2328 pvar = pfi->pvar1;
2329
2330 /* Grid data */
2331 if (pfi->type == 1)
2332 {
2333 pfi->gsiz = pfi->dnum[0] * pfi->dnum[1];
2334 /* if (pfi->ppflag) pfi->gsiz = pfi->ppisiz * pfi->ppjsiz; */
2335 /* add the XY header to gsiz */
2336 if (pfi->xyhdr)
2337 {
2338 if (pvar->dfrm == 1)
2339 {
2340 pfi->xyhdr = pfi->xyhdr * 4 / 1;
2341 }
2342 else if (pvar->dfrm == 2 || pvar->dfrm == -2)
2343 {
2344 pfi->xyhdr = pfi->xyhdr * 4 / 2;
2345 }
2346 else if (pfi->flt64)
2347 {
2348 pfi->xyhdr = pfi->xyhdr * 4 / 8;
2349 }
2350 pfi->gsiz = pfi->gsiz + pfi->xyhdr;
2351 }
2352
2353 /* adjust the size of hdrb and trlb for non-float data */
2354 if (pvar->dfrm == 1)
2355 {
2356 hdrb = hdrb * 4 / 1;
2357 trlb = trlb * 4 / 1;
2358 }
2359 else if (pvar->dfrm == 2 || pvar->dfrm == -2)
2360 {
2361 hdrb = hdrb * 4 / 2;
2362 trlb = trlb * 4 / 2;
2363 }
2364
2365 if (pfi->seqflg)
2366 {
2367 /* pad the grid size with 2 4-byte chunks */
2368 if (pvar->dfrm == 1)
2369 {
2370 pfi->gsiz += 8;
2371 }
2372 else if (pvar->dfrm == 2 || pvar->dfrm == -2)
2373 {
2374 pfi->gsiz += 4;
2375 }
2376 else
2377 {
2378 if (pfi->flt64)
2379 pfi->gsiz += 1;
2380 else
2381 pfi->gsiz += 2;
2382 }
2383 /* pad the header with 2 4-byte chunks*/
2384 if (hdrb > 0)
2385 {
2386 if (pvar->dfrm == 1)
2387 {
2388 hdrb = hdrb + 8;
2389 }
2390 else if (pvar->dfrm == 2 || pvar->dfrm == -2)
2391 {
2392 hdrb = hdrb + 4;
2393 }
2394 else
2395 {
2396 hdrb += 2;
2397 }
2398 }
2399 /* how far we have to go into the file before getting to 1st var */
2400 if (pvar->dfrm == 1)
2401 {
2402 pvar->offset = 4 + hdrb;
2403 acum = 4 + hdrb;
2404 }
2405 else if (pvar->dfrm == 2 || pvar->dfrm == -2)
2406 {
2407 pvar->offset = 2 + hdrb;
2408 acum = 2 + hdrb;
2409 }
2410 else
2411 {
2412 pvar->offset = 1 + hdrb;
2413 acum = 1 + hdrb;
2414 }
2415 }
2416 else
2417 {
2418 /* how far we have to go into the file before getting to 1st var */
2419 pvar->offset = hdrb;
2420 acum = hdrb;
2421 }
2422
2423 levs = pvar->levels;
2424 if (levs == 0) levs = 1;
2425 pvar->recoff = 0;
2426 recacm = 0;
2427 pvar++;
2428
2429 for (i = 1; i < pfi->vnum; i++)
2430 {
2431 if (pvar->var_t)
2432 {
2433 acum = acum + levs * (pfi->gsiz) * (pfi->dnum[3]);
2434 }
2435 else
2436 {
2437 acum = acum + (levs * pfi->gsiz);
2438 // acumstride = acum ;
2439 }
2440 recacm += levs;
2441 pvar->offset = acum;
2442 pvar->recoff = recacm;
2443 levs = pvar->levels;
2444 if (levs == 0) levs = 1;
2445 pvar++;
2446 }
2447
2448 recacm += levs;
2449
2450 /* last variable */
2451 acum = acum + (levs * pfi->gsiz);
2452
2453 pfi->tsiz = acum;
2454 pfi->trecs = recacm;
2455 if (pfi->seqflg) pfi->tsiz -= 1;
2456 pfi->tsiz += trlb;
2457 }
2458 else
2459 {
2460 fprintf(stderr, "Grid data type unsupported!");
2461 return -1;
2462 }
2463
2464 /* set the global calendar and check if we are trying to change with a new
2465 file... we do this here to set the calandar for templating */
2466 if (/*mfcmn.*/ cal365 < 0)
2467 {
2468 /*mfcmn.*/ cal365 = pfi->calendar;
2469 }
2470 else
2471 {
2472 if (pfi->calendar != /*mfcmn.*/ cal365)
2473 {
2474 gaprnt(0, "Attempt to change the global calendar...\n");
2475 if (/*mfcmn.*/ cal365)
2476 {
2477 gaprnt(0, "The calendar is NOW 365 DAYS and you attempted to "
2478 "open a standard calendar file\n");
2479 }
2480 else
2481 {
2482 gaprnt(0, "The calendar is NOW STANDARD and you attempted to "
2483 "open a 365-day calendar file\n");
2484 }
2485 goto retrn;
2486 }
2487 }
2488
2489 /* If the file name is a time series template, figure out
2490 which times go with which files, so we don't waste a lot
2491 of time later opening and closing files unnecessarily. */
2492
2493 if (pfi->tmplat)
2494 {
2495 /* The fnums array is the size of the time axis
2496 multiplied by the size of the ensemble axis.
2497 It contains the t index which generates the filename
2498 that contains the data for each timestep.
2499 If the ensemble has no data file for a given time,
2500 the fnums value will be -1 */
2501 pfi->fnums = (gaint *) galloc(sizeof(gaint) * pfi->dnum[3] * pfi->dnum[4], "fnums1");
2502 if (pfi->fnums == NULL)
2503 {
2504 gaprnt(0, "Open Error: memory allocation failed for fnums\n");
2505 goto err8;
2506 }
2507 /* get dt structure for t=1 */
2508 gr2t(pfi->grvals[3], 1.0, &tdefi);
2509 /* loop over ensembles */
2510 ens = pfi->ens1;
2511 e = 1;
2512 while (e <= pfi->dnum[4])
2513 {
2514 j = -1;
2515 t = 1;
2516 /* set fnums value to -1 for time steps before ensemble initial time
2517 */
2518 while (t < ens->gt)
2519 {
2520 pfi->fnums[t - 1] = j;
2521 t++;
2522 }
2523 j = ens->gt;
2524 /* get dt structure for ensemble initial time */
2525 gr2t(pfi->grvals[3], ens->gt, &tdefe);
2526 /* get filename for initial time of current ensemble member */
2527 ch = gafndt(pfi->name, &tdefe, &tdefe, pfi->abvals[3], pfi->pchsub1, pfi->ens1, ens->gt, e, &flag);
2528 if (ch == NULL)
2529 {
2530 sprintf(pout, "Open Error: couldn't determine data file name for e=%d t=%d\n", e, ens->gt);
2531 gaprnt(0, pout);
2532 goto err8;
2533 }
2534 /* set the pfi->tmplat flag to the flag returned by gafndt */
2535 if (flag == 0)
2536 {
2537 gaprnt(1, "Warning: OPTIONS keyword \"template\" is used, but the \n");
2538 gaprnt(1, " DSET entry contains no substitution templates.\n");
2539 pfi->tmplat = 1;
2540 }
2541 else
2542 {
2543 pfi->tmplat = flag;
2544 }
2545 /* for non-indexed, non-netcdf/hdf, gridded data */
2546 if (pfi->type == 1)
2547 { /* gridded data */
2548 if (pfi->ncflg == 0)
2549 { /* not netcdf/hdf */
2550 if (pfi->idxflg == 0)
2551 { /* not indexed */
2552 if ((flag == 1) && (pfi->dnum[4] > 1))
2553 {
2554 gaprnt(0, "Open Error: If the data type is gridded binary, \n");
2555 gaprnt(0, " and the E dimension size is greater than 1 \n");
2556 gaprnt(0, " and templating in the T dimension is used,\n");
2557 gaprnt(0, " then templating in the E dimension must also be used.\n");
2558 goto retrn;
2559 }
2560 }
2561 else if (pfi->idxflg == 1)
2562 { /* GRIB1 */
2563 if ((flag < 2) && (pfi->dnum[4] > 1))
2564 {
2565 gaprnt(0, "Open Error: If the data type is GRIB1 \n");
2566 gaprnt(0, " and the E dimension size is greater than 1 \n");
2567 gaprnt(0, " then templating in the E dimension must be used.\n");
2568 goto retrn;
2569 }
2570 }
2571 }
2572 }
2573 pfi->fnums[t - 1] = j;
2574 /* loop over remaining valid times for this ensemble */
2575 for (t = ens->gt + 1; t < ens->gt + ens->length; t++)
2576 {
2577 /* get filename for time index=t ens=e */
2578 gr2t(pfi->grvals[3], (gadouble) t, &tdef);
2579 pos = gafndt(pfi->name, &tdef, &tdefe, pfi->abvals[3], pfi->pchsub1, pfi->ens1, t, e, &flag);
2580 if (pos == NULL)
2581 {
2582 sprintf(pout, "Open Error: couldn't determine data file name for e=%d t=%d\n", e, t);
2583 gaprnt(0, pout);
2584 goto err8;
2585 }
2586 if (strcmp(ch, pos) != 0)
2587 { /* filename has changed */
2588 j = t;
2589 gree(ch, "f176");
2590 ch = pos;
2591 }
2592 else
2593 {
2594 gree(pos, "f176a");
2595 }
2596 pfi->fnums[+t - 1] = j;
2597 }
2598 gree(ch, "f177");
2599
2600 /* set fnums value to -1 for time steps after ensemble final time */
2601 j = -1;
2602 while (t <= pfi->dnum[3])
2603 {
2604 pfi->fnums[t - 1] = j;
2605 t++;
2606 }
2607 e++;
2608 ens++;
2609 }
2610 pfi->fnumc = 0;
2611 pfi->fnume = 0;
2612 }
2613
2614 fclose(descr);
2615
2616 return (status);
2617
2618 err1:
2619 gaprnt(0, "Open Error: Missing or invalid dimension size.\n");
2620 goto err9;
2621
2622 err2:
2623 gaprnt(0, "Open Error: Missing or invalid dimension");
2624 gaprnt(0, " scaling type\n");
2625 goto err9;
2626
2627 err3a_tdef:
2628 gaprnt(0, "Open Error: Start Time missing in tdef card");
2629 gaprnt(0, " starting value\n");
2630 goto err9;
2631
2632 err3b_tdef:
2633 gaprnt(0, "Open Error: Invalid start time in tdef card");
2634 gaprnt(0, " starting value\n");
2635 goto err9;
2636
2637 err3c_tdef:
2638 gaprnt(0, "Open Error: Missing or invalid dimension");
2639 gaprnt(0, " starting value\n");
2640 goto err9;
2641
2642 err4a_tdef:
2643 gaprnt(0, "Open Error: Time increment missing in tdef\n");
2644 gaprnt(0, " use 1 for single time data\n");
2645 goto err9;
2646
2647 err4b_tdef:
2648 gaprnt(0, "Open Error: Invalid time increment in tdef\n");
2649 gaprnt(0, " use 1 for single time data\n");
2650 goto err9;
2651
2652 err4c_tdef:
2653 gaprnt(0, "Open Error: 0 time increment in tdef\n");
2654 gaprnt(0, " use 1 for single time data\n");
2655 goto err9;
2656
2657 err5:
2658 gaprnt(0, "Open Error: Missing or invalid variable");
2659 gaprnt(0, " count\n");
2660 goto err9;
2661
2662 err6:
2663 gaprnt(0, "Open Error: Invalid variable record\n");
2664 goto err9;
2665
2666 err8:
2667 gaprnt(0, "Open Error: Memory allocation Error in gaddes.c\n");
2668 goto retrn;
2669
2670 err9:
2671 gaprnt(0, " --> The invalid description file record is: \n");
2672 gaprnt(0, " --> ");
2673 gaprnt(0, mrec);
2674 gaprnt(0, "\n");
2675
2676 retrn:
2677 gaprnt(0, " The data file was not opened. \n");
2678 fclose(descr);
2679 return (1);
2680 }
2681