1 /* Copyright (C) 1988-2005 by Brian Doty and the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 /* Authored by B. Doty */
7
8 /* kk 020624 --- Change for 64bit seek K.Komine */
9
10 #ifdef HAVE_CONFIG_H
11 #include <config.h>
12
13 /* If autoconfed, only include malloc.h when it's present */
14 #ifdef HAVE_MALLOC_H
15 #include <malloc.h>
16 #endif
17
18 #else /* undef HAVE_CONFIG_H */
19
20 #include <malloc.h>
21
22 #endif /* HAVE_CONFIG_H */
23
24
25 #include <stdio.h>
26 #include <string.h>
27 #include <math.h>
28 #include <ctype.h>
29 #include "grads.h"
30
31 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
32 extern struct gamfcmn mfcmn;
33 /*mf 961205 --- expose Mike Fiorino's global struct to these routines for 365 day calandars mf*/
34
35 static char pout[512];
36 static FILE *pdfi; /* File descriptor for pdef file */
37
38 FILE *descr; /* File descriptor pointer */
39 /* Also used, via extern, by gasdf.c */
40
41 void ll2eg (int , int, float *, float, float, float *, float *, float *) ;
42 void ll2pse(int, int, float *, float, float, float *, float *);
43 void ll2ops (float *, float , float , float *, float *) ;
44
45
46 /* Read a GrADS data descriptor file and fill in the information
47 into a gafile structure. The gafile structure should be
48 allocated and must be initialized by getpfi. If this routine
49 returns an error, release the pfi structure and allocated
50 storage via frepfi. mflag indicates whether to read the
51 stnmap/index file; if this routine is being called to
52 preprocess the dd file then this flag should be 0.
53 A mflag value of 2 will check if the stnmap/index file can be
54 opened, but will not read it; mflag==2 will also turn off
55 calculation of the pdef interpolation values */
56
gaddes(char * name,struct gafile * pfi,int mflag)57 int gaddes (char *name, struct gafile *pfi, int mflag) {
58 struct gavar *pvar,*pvar2;
59 struct dt tdef,tdefi,dt1,dt2;
60 struct gaindx *pindx;
61 struct gaattr *attrib, *testattr;
62 struct gachsub *pchsub;
63 float *vals;
64 int size=0,rc,len,swpflg,cnt,flag,tim1,tim2;
65 char rec[512], mrec[512], *ch, *pos, *sname, *vectorpairs, *pair;
66 int flgs[8],cflg,i,j,ii,jj,err,hdrb,trlb,mflflg;
67 int mcnt,maxlv,foundvar1,foundvar2;
68 off_t levs,acum,acumvz,recacm;
69 float temp,v1,v2;
70 FILE *mfile;
71 char pdefnm[256],var1[256],var2[256];
72 char *varname,*attrname,*attrtype,attrvalue[512];
73 int pdefop1=0, pdefop2=0;
74 int nzstride=0,first=1;
75 int levsua=0,acumstride=0,npairs;
76 unsigned char vermap;
77 int idum,reclen;
78 float fdum;
79 unsigned char urec[4];
80 int diag=0;
81 #if GRADS_CRAY == 1
82 int crayflg=0;
83 #endif
84
85 static char *errs[9] = {"XDEF","YDEF","ZDEF","TDEF","UNDEF",
86 "DSET","VARS","TITLE","DTYPE"};
87
88 /*mf --- define here vice grads.c for cdunif.c mf*/
89 mfcmn.fullyear=-999;
90
91 /* initialize variables */
92 hdrb = 0;
93 trlb = 0;
94 pdfi = NULL;
95 mflflg = 0;
96 mfile = NULL;
97 pfi->mfile = NULL;
98 mcnt = -1;
99 vectorpairs = NULL;
100 pair = NULL;
101 varname = attrname = attrtype = NULL;
102 attrib = NULL;
103 sname=NULL;
104
105 /* Try to open descriptor file */
106 descr = fopen (name, "r");
107 if (descr == NULL) {
108 /* Try adding default suffix of .ctl */
109 sname = (char *)malloc(strlen(name)+5);
110 if(sname == NULL) {
111 gaprnt(0,"malloc error in creating date descriptor file name\n");
112 return(1);
113 }
114 for(i=0;i<=strlen(name);i++) *(sname+i)=*(name+i);
115 strcat(sname,".ctl");
116 descr = fopen (sname, "r");
117 }
118
119 /* If still can't open descriptor file, give up */
120 if (descr == NULL) {
121 gaprnt (0,"Open Error: Can't open description file\n");
122 if (sname) free (sname);
123 return(1);
124 }
125
126 /* Copy descriptor file name into gafile structure */
127 if (sname != NULL) {
128 getwrd (pfi->dnam,sname,512);
129 free(sname);
130 } else {
131 getwrd (pfi->dnam,name,512);
132 }
133
134 /* initialize error flags */
135 for (i=0;i<8;i++) flgs[i] = 1;
136
137 /* Parse the data descriptor file */
138 while (fgets(rec,512,descr)!=NULL) {
139
140 /* Remove any leading blanks from rec */
141 reclen = strlen(rec);
142 jj = 0;
143 while (jj<reclen && rec[0]==' ') {
144 for (ii=0; ii<reclen; ii++) rec[ii] = rec[ii+1];
145 jj++;
146 }
147
148 /* Keep mixed case and lower case versions of rec handy */
149 strcpy (mrec,rec);
150 lowcas(rec);
151
152 if ( (cmpwrd("*",rec)) || (cmpwrd("#",rec)) || !isalnum(rec[0]) ) {
153 cflg = 1;
154 /* Parse comment if it contains attribute metadata */
155 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0)) {
156 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
157 }
158
159 } else if (cmpwrd("byteswapped",rec)) {
160 pfi->bswap = 1;
161
162 } else if (cmpwrd("fileheader",rec)) {
163 if ( (ch=nxtwrd(rec))==NULL ) {
164 gaprnt (1,"Description file warning: Missing fileheader length\n");
165 } else {
166 ch = longprs(ch,&(pfi->fhdr));
167 if (ch==NULL) {
168 gaprnt (1,"Fileheader record invalid\n");
169 pfi->fhdr = 0;
170 }
171 }
172
173 } else if (cmpwrd("xyheader",rec)) {
174 if ( (ch=nxtwrd(rec))==NULL ) {
175 gaprnt (1,"Description file warning: Missing xy grid header length\n");
176 } else {
177 ch = longprs(ch,&(pfi->xyhdr));
178 if (ch==NULL) {
179 gaprnt (1,"xy grid header length invalid\n");
180 pfi->xyhdr = 0;
181 } else {
182 pfi->xyhdr = pfi->xyhdr / sizeof(float);
183 }
184 }
185
186 } else if (cmpwrd("unpack",rec)) {
187 if ( (ch=nxtwrd(mrec))==NULL ) {
188 gaprnt (1,"Descriptor File Warning: Missing attribute names in unpack record\n");
189 } else {
190 /* get the scale factor attribute name */
191 len = 0;
192 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
193 pfi->scattr = (char *)malloc(len+3);
194 for (i=0; i<len; i++) *(pfi->scattr+i) = *(ch+i);
195 *(pfi->scattr+len) = '\0';
196 /* set the packflg to 1, meaning only scale factor has been retrieved */
197 pfi->packflg = 1;
198
199 /* get the offset attribute name */
200 if ( (ch=nxtwrd(ch)) == NULL ) {
201 gaprnt (2,"Descriptor File Warning: No offset attribute name in unpack record\n");
202 } else {
203 len = 0;
204 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
205 pfi->ofattr = (char *)malloc(len+3);
206 for (i=0; i<len; i++) *(pfi->ofattr+i) = *(ch+i);
207 *(pfi->ofattr+len) = '\0';
208 /* Set the packflg to 2, meaning scale factor and offset have been retrieved */
209 pfi->packflg = 2;
210 }
211 }
212
213 } else if (cmpwrd("theader",rec)) {
214 if ( (ch=nxtwrd(rec))==NULL ) {
215 gaprnt (1,"Description file warning: Missing time chunk header length\n");
216 } else {
217 ch = longprs(ch,&(pfi->thdr));
218 if (ch==NULL) {
219 gaprnt (1,"time chunk grid header length invalid\n");
220 pfi->thdr = 0;
221 } else {
222 pfi->thdr = pfi->thdr / sizeof(float);
223 }
224 }
225
226 } else if (cmpwrd("format",rec) || cmpwrd("options",rec)) {
227 if ( (ch=nxtwrd(rec))==NULL ) {
228 gaprnt (1,"Description file warning: Missing options keyword\n");
229 } else {
230 while (ch!=NULL) {
231 if (cmpwrd("sequential",ch)) pfi->seqflg = 1;
232 else if (cmpwrd("yrev",ch)) pfi->yrflg = 1;
233 else if (cmpwrd("zrev",ch)) pfi->zrflg = 1;
234 else if (cmpwrd("template",ch)) pfi->tmplat = 1;
235 else if (cmpwrd("byteswapped",ch)) pfi->bswap = 1;
236 else if (cmpwrd("cray_32bit_ieee",ch) && GRADS_CRAY) pfi->cray_ieee = 1;
237 else if (cmpwrd("cray_32bit_ieee",ch) && !GRADS_CRAY) pfi->cray_ieee = 0;
238 else if (cmpwrd("365_day_calendar",ch)) {
239 pfi->calendar=1;
240 mfcmn.cal365=pfi->calendar;
241 }
242 else if (cmpwrd("big_endian",ch)) {
243 if (!BYTEORDER) pfi->bswap = 1;
244 }
245 else if (cmpwrd("little_endian",ch)) {
246 if (BYTEORDER) pfi->bswap = 1;
247 }
248 else {
249 gaprnt (0,"Open Error: Data file type invalid\n");
250 goto err9;
251 }
252 ch = nxtwrd(ch);
253 }
254 }
255
256 } else if (cmpwrd("trailerbytes",rec)) {
257 if ( (ch=nxtwrd(rec))==NULL ) {
258 gaprnt (1,"Trailerbytes record invalid\n");
259 } else {
260 ch = intprs(ch,&trlb);
261 if (ch==NULL) {
262 gaprnt (1,"Trailerbytes record invalid\n");
263 trlb = 0;
264 } else {
265 trlb = trlb/4;
266 }
267 }
268
269
270 } else if (cmpwrd("headerbytes",rec)) {
271 if ( (ch=nxtwrd(rec))==NULL ) {
272 gaprnt (1,"Headerbytes record invalid\n");
273 } else {
274 ch = intprs(ch,&hdrb);
275 if (ch==NULL) {
276 gaprnt (1,"Headerbytes record invalid\n");
277 hdrb = 0;
278 } else {
279 hdrb = hdrb/4;
280 }
281 }
282
283 /* Handle the chsub records. time1, time2, then a string, multiple times */
284 } else if (cmpwrd("chsub",rec)) {
285
286 pchsub = pfi->pchsub1; /* point to first block in chain */
287 if (pchsub) {
288 while (pchsub->forw!=NULL) {
289 pchsub = pchsub->forw; /* advance to end of chain */
290 }
291 }
292 flag = 0;
293 ch = mrec;
294 while (1) {
295 if ( (ch=nxtwrd(ch)) == NULL ) break;
296 flag = 1;
297 if ( (ch = intprs(ch,&tim1)) == NULL) break;
298 if ( (ch=nxtwrd(ch)) == NULL ) break;
299 if (*ch=='*' && (*(ch+1)==' '||*(ch+1)=='\t')) tim2 = -99;
300 else if ( (ch = intprs(ch,&tim2)) == NULL) break;
301 if ( (ch=nxtwrd(ch)) == NULL ) break;
302 flag = 0;
303 if (pchsub) {
304 pchsub->forw = (struct gachsub *)calloc(1,sizeof(struct gachsub));
305 if (pchsub->forw==NULL) goto err8;
306 pchsub = pchsub->forw;
307 } else {
308 pfi->pchsub1 = (struct gachsub *)calloc(1,sizeof(struct gachsub));
309 if (pfi->pchsub1==NULL) goto err8;
310 pchsub = pfi->pchsub1;
311 }
312 len = wrdlen(ch);
313 pchsub->ch = (char *)malloc(len+1);
314 if (pchsub->ch==NULL) goto err8;
315 getwrd(pchsub->ch,ch,len);
316 pchsub->t1 = tim1;
317 pchsub->t2 = tim2;
318 }
319 if (flag) {
320 gaprnt (1,"Description file warning: Invalid chsub record; Ignored\n");
321 }
322
323 } else if (cmpwrd("dtype",rec)) {
324 if ( (ch=nxtwrd(rec))==NULL ) {
325 gaprnt (1,"Description file warning: Missing data type\n");
326 gaprnt (1," Assuming data file type is grid\n");
327 }
328 if (cmpwrd("station",ch)) {
329 pfi->type = 2;
330 flgs[0] = 0;
331 flgs[1] = 0;
332 flgs[2] = 0;
333
334 } else if (cmpwrd("bufr",ch)) {
335 pfi->type = 2; /* station data type */
336 mflag = 0; /* don't try to read a stnmap file */
337 pfi->bufrflg = 1;
338 flgs[0] = 0;
339 flgs[1] = 0;
340 flgs[2] = 0;
341 /* allocate memory for the bufrinfo structure and the two bufrtimeinfo structures */
342 pfi->bufrinfo = (struct bufrinfo *)malloc(sizeof(struct bufrinfo));
343 /* initialize with bad values */
344 for (j=0;j<2;j++) {
345 pfi->bufrinfo->lonxy[j] = pfi->bufrinfo->latxy[j] = -999;;
346 pfi->bufrinfo->levxy[j] = pfi->bufrinfo->stidxy[j] = -999;;
347 pfi->bufrinfo->base.yrxy[j] = pfi->bufrinfo->base.moxy[j] = -999;;
348 pfi->bufrinfo->base.dyxy[j] = pfi->bufrinfo->base.hrxy[j] = -999;;
349 pfi->bufrinfo->base.mnxy[j] = pfi->bufrinfo->offset.yrxy[j] = -999;;
350 pfi->bufrinfo->offset.moxy[j] = pfi->bufrinfo->offset.dyxy[j] = -999;;
351 pfi->bufrinfo->offset.hrxy[j] = pfi->bufrinfo->offset.mnxy[j] = -999;;
352 }
353
354 } else if (cmpwrd("grib",ch)) {
355 pfi->idxflg = 1;
356 if ( (ch=nxtwrd(ch))!=NULL ) {
357 if ( intprs(ch,&(pfi->grbgrd))==NULL) {
358 gaprnt (1,"Description file warning: Invalid GRIB option\n");
359 pfi->grbgrd = -999;
360 }
361 }
362 }
363 #if USESDF ==1
364 else if (cmpwrd("netcdf",ch)) pfi->ncflg = 1;
365 #if USEHDF ==1
366 else if (cmpwrd("hdfsds",ch)) pfi->ncflg = 2;
367 #endif
368 #endif
369 else {
370 gaprnt (0,"Open Error: Data file type invalid\n");
371 goto err9;
372 }
373
374 } else if (cmpwrd("xvar",rec)) {
375 if ( (ch=nxtwrd(rec))==NULL ) {
376 gaprnt (1,"Description file warning: Missing x,y pair for XVAR entry\n");
377 } else {
378 j = 0;
379 while (1) {
380 if ( (ch=intprs(ch,&(pfi->bufrinfo->lonxy[j])))==NULL ) goto err6a;
381 while (*ch==' ') ch++;
382 if (*ch!=',') break;
383 ch++;
384 while (*ch==' ') ch++;
385 j++;
386 if (j>1) goto err6a;
387 }
388 if (pfi->bufrinfo->lonxy[0]==-999 || pfi->bufrinfo->lonxy[1]==-999) goto err6a;
389 }
390
391 } else if (cmpwrd("yvar",rec)) {
392 if ( (ch=nxtwrd(rec))==NULL ) {
393 gaprnt (1,"Description file warning: Missing x,y pair for YVAR entry\n");
394 } else {
395 j = 0;
396 while (1) {
397 if ( (ch=intprs(ch,&(pfi->bufrinfo->latxy[j])))==NULL ) goto err6a;
398 while (*ch==' ') ch++;
399 if (*ch!=',') break;
400 ch++;
401 while (*ch==' ') ch++;
402 j++;
403 if (j>1) goto err6a;
404 }
405 if (pfi->bufrinfo->latxy[0]==-999 || pfi->bufrinfo->latxy[1]==-999) goto err6a;
406 }
407
408 } else if (cmpwrd("zvar",rec)) {
409 if ( (ch=nxtwrd(rec))==NULL ) {
410 gaprnt (1,"Description file warning: Missing x,y pair for ZVAR entry\n");
411 } else {
412 j = 0;
413 while (1) {
414 if ( (ch=intprs(ch,&(pfi->bufrinfo->levxy[j])))==NULL ) goto err6a;
415 while (*ch==' ') ch++;
416 if (*ch!=',') break;
417 ch++;
418 while (*ch==' ') ch++;
419 j++;
420 if (j>1) goto err6a;
421 }
422 if (pfi->bufrinfo->levxy[0]==-999 || pfi->bufrinfo->levxy[1]==-999) goto err6a;
423 }
424
425 } else if (cmpwrd("tvar",rec)) {
426 if ( (ch=nxtwrd(rec))==NULL ) {
427 gaprnt (1,"Description file warning: Missing x,y pairs for TVAR entry\n");
428 } else {
429 while (nxtwrd(ch)!=NULL) {
430 if (cmpwrd("yr",ch)) {
431 if ((ch=nxtwrd(ch))==NULL) {
432 gaprnt (1,"Description file warning: Missing x,y pair for TVAR yr entry\n");
433 } else {
434 j = 0;
435 while (1) {
436 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.yrxy[j])))==NULL ) goto err6a;
437 while (*ch==' ') ch++;
438 if (*ch!=',') break;
439 ch++;
440 while (*ch==' ') ch++;
441 j++;
442 if (j>1) goto err6a;
443 }
444 if (pfi->bufrinfo->base.yrxy[0]==-999 || pfi->bufrinfo->base.yrxy[1]==-999) goto err6a;
445 }
446
447 } else if (cmpwrd("mo",ch)) {
448 if ((ch=nxtwrd(ch))==NULL) {
449 gaprnt (1,"Description file warning: Missing x,y pair for TVAR mo entry\n");
450 } else {
451 j = 0;
452 while (1) {
453 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.moxy[j])))==NULL ) goto err6a;
454 while (*ch==' ') ch++;
455 if (*ch!=',') break;
456 ch++;
457 while (*ch==' ') ch++;
458 j++;
459 if (j>1) goto err6a;
460 }
461 if (pfi->bufrinfo->base.moxy[0]==-999 || pfi->bufrinfo->base.moxy[1]==-999) goto err6a;
462 }
463 } else if (cmpwrd("dy",ch)) {
464 if ((ch=nxtwrd(ch))==NULL) {
465 gaprnt (1,"Description file warning: Missing x,y pair for TVAR dy entry\n");
466 } else {
467 j = 0;
468 while (1) {
469 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.dyxy[j])))==NULL ) goto err6a;
470 while (*ch==' ') ch++;
471 if (*ch!=',') break;
472 ch++;
473 while (*ch==' ') ch++;
474 j++;
475 if (j>1) goto err6a;
476 }
477 if (pfi->bufrinfo->base.dyxy[0]==-999 || pfi->bufrinfo->base.dyxy[1]==-999) goto err6a;
478 }
479 } else if (cmpwrd("hr",ch)) {
480 if ((ch=nxtwrd(ch))==NULL) {
481 gaprnt (1,"Description file warning: Missing x,y pair for TVAR hr entry\n");
482 } else {
483 j = 0;
484 while (1) {
485 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.hrxy[j])))==NULL ) goto err6a;
486 while (*ch==' ') ch++;
487 if (*ch!=',') break;
488 ch++;
489 while (*ch==' ') ch++;
490 j++;
491 if (j>1) goto err6a;
492 }
493 if (pfi->bufrinfo->base.hrxy[0]==-999 || pfi->bufrinfo->base.hrxy[1]==-999) goto err6a;
494 }
495 } else if (cmpwrd("mn",ch)) {
496 if ((ch=nxtwrd(ch))==NULL) {
497 gaprnt (1,"Description file warning: Missing x,y pair for TVAR mn entry\n");
498 } else {
499 j = 0;
500 while (1) {
501 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.mnxy[j])))==NULL ) goto err6a;
502 while (*ch==' ') ch++;
503 if (*ch!=',') break;
504 ch++;
505 while (*ch==' ') ch++;
506 j++;
507 if (j>1) goto err6a;
508 }
509 if (pfi->bufrinfo->base.mnxy[0]==-999 || pfi->bufrinfo->base.mnxy[1]==-999) goto err6a;
510 }
511 } else if (cmpwrd("sc",ch)) {
512 if ((ch=nxtwrd(ch))==NULL) {
513 gaprnt (1,"Description file warning: Missing x,y pair for TVAR sc entry\n");
514 } else {
515 j = 0;
516 while (1) {
517 if ( (ch=intprs(ch,&(pfi->bufrinfo->base.scxy[j])))==NULL ) goto err6a;
518 while (*ch==' ') ch++;
519 if (*ch!=',') break;
520 ch++;
521 while (*ch==' ') ch++;
522 j++;
523 if (j>1) goto err6a;
524 }
525 if (pfi->bufrinfo->base.scxy[0]==-999 || pfi->bufrinfo->base.scxy[1]==-999) goto err6a;
526 }
527 } else {
528 goto err6a;
529 }
530 } /* end of while loop */
531 }
532
533 } else if (cmpwrd("toffvar",rec)) {
534 if ( (ch=nxtwrd(rec))==NULL ) {
535 gaprnt (1,"Description file warning: Missing x,y pairs for TOFFVAR entry\n");
536 } else {
537 while (nxtwrd(ch)!=NULL) {
538 if (cmpwrd("yr",ch)) {
539 if ((ch=nxtwrd(ch))==NULL) {
540 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR yr entry\n");
541 } else {
542 j = 0;
543 while (1) {
544 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.yrxy[j])))==NULL ) goto err6a;
545 while (*ch==' ') ch++;
546 if (*ch!=',') break;
547 ch++;
548 while (*ch==' ') ch++;
549 j++;
550 if (j>1) goto err6a;
551 }
552 if (pfi->bufrinfo->offset.yrxy[0]==-999 || pfi->bufrinfo->offset.yrxy[1]==-999) goto err6a;
553 }
554
555 } else if (cmpwrd("mo",ch)) {
556 if ((ch=nxtwrd(ch))==NULL) {
557 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR mo entry\n");
558 } else {
559 j = 0;
560 while (1) {
561 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.moxy[j])))==NULL ) goto err6a;
562 while (*ch==' ') ch++;
563 if (*ch!=',') break;
564 ch++;
565 while (*ch==' ') ch++;
566 j++;
567 if (j>1) goto err6a;
568 }
569 if (pfi->bufrinfo->offset.moxy[0]==-999 || pfi->bufrinfo->offset.moxy[1]==-999) goto err6a;
570 }
571 } else if (cmpwrd("dy",ch)) {
572 if ((ch=nxtwrd(ch))==NULL) {
573 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR dy entry\n");
574 } else {
575 j = 0;
576 while (1) {
577 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.dyxy[j])))==NULL ) goto err6a;
578 while (*ch==' ') ch++;
579 if (*ch!=',') break;
580 ch++;
581 while (*ch==' ') ch++;
582 j++;
583 if (j>1) goto err6a;
584 }
585 if (pfi->bufrinfo->offset.dyxy[0]==-999 || pfi->bufrinfo->offset.dyxy[1]==-999) goto err6a;
586 }
587 } else if (cmpwrd("hr",ch)) {
588 if ((ch=nxtwrd(ch))==NULL) {
589 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR hr entry\n");
590 } else {
591 j = 0;
592 while (1) {
593 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.hrxy[j])))==NULL ) goto err6a;
594 while (*ch==' ') ch++;
595 if (*ch!=',') break;
596 ch++;
597 while (*ch==' ') ch++;
598 j++;
599 if (j>1) goto err6a;
600 }
601 if (pfi->bufrinfo->offset.hrxy[0]==-999 || pfi->bufrinfo->offset.hrxy[1]==-999) goto err6a;
602 }
603 } else if (cmpwrd("mn",ch)) {
604 if ((ch=nxtwrd(ch))==NULL) {
605 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR mn entry\n");
606 } else {
607 j = 0;
608 while (1) {
609 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.mnxy[j])))==NULL ) goto err6a;
610 while (*ch==' ') ch++;
611 if (*ch!=',') break;
612 ch++;
613 while (*ch==' ') ch++;
614 j++;
615 if (j>1) goto err6a;
616 }
617 if (pfi->bufrinfo->offset.mnxy[0]==-999 || pfi->bufrinfo->offset.mnxy[1]==-999) goto err6a;
618 }
619 } else if (cmpwrd("sc",ch)) {
620 if ((ch=nxtwrd(ch))==NULL) {
621 gaprnt (1,"Description file warning: Missing x,y pair for TOFFVAR sc entry\n");
622 } else {
623 j = 0;
624 while (1) {
625 if ( (ch=intprs(ch,&(pfi->bufrinfo->offset.scxy[j])))==NULL ) goto err6a;
626 while (*ch==' ') ch++;
627 if (*ch!=',') break;
628 ch++;
629 while (*ch==' ') ch++;
630 j++;
631 if (j>1) goto err6a;
632 }
633 if (pfi->bufrinfo->offset.scxy[0]==-999 || pfi->bufrinfo->offset.scxy[1]==-999) goto err6a;
634 }
635 } else {
636 goto err6a;
637 }
638 } /* end of while loop */
639 }
640
641 } else if (cmpwrd("stid",rec)) {
642 if ( (ch=nxtwrd(rec))==NULL ) {
643 gaprnt (1,"Description file warning: Missing x,y pair for STID entry\n");
644 } else {
645 j = 0;
646 while (1) {
647 if ( (ch=intprs(ch,&(pfi->bufrinfo->stidxy[j])))==NULL ) goto err6a;
648 while (*ch==' ') ch++;
649 if (*ch!=',') break;
650 ch++;
651 while (*ch==' ') ch++;
652 j++;
653 if (j>1) goto err6a;
654 }
655 if (pfi->bufrinfo->stidxy[0]==-999 || pfi->bufrinfo->stidxy[1]==-999) goto err6a;
656 }
657
658 } else if (cmpwrd("title",rec)) {
659 if ( (ch=nxtwrd(mrec))==NULL ) {
660 gaprnt (1,"Description file warning: Missing title string\n");
661 } else {
662 getstr (pfi->title,ch,512);
663 flgs[7] = 0;
664 }
665
666 } else if (cmpwrd("dset",rec)) {
667 ch = nxtwrd(mrec);
668 if (ch==NULL) {
669 gaprnt (0,"Descriptor File Error: Data file name is missing\n");
670 goto err9;
671 }
672 if (*ch=='^' || *ch=='$') {
673 fnmexp (pfi->name,ch,name);
674 } else {
675 getwrd (pfi->name,ch,512);
676 }
677 flgs[5] = 0;
678
679 } else if (cmpwrd("stnmap",rec) || cmpwrd("index",rec)) {
680 /* map file processing */
681 if (cmpwrd("index",rec)) pfi->idxflg = 1;
682 ch = nxtwrd(mrec);
683 if (ch==NULL) {
684 gaprnt (0,"Open Error: Station or Index Map file name is missing\n");
685 goto err9;
686 }
687 if (*ch=='^' || *ch=='$') {
688 fnmexp (pout, ch, name);
689 } else {
690 getwrd (pout, ch, 500);
691 }
692 len = 0;
693 while (*(pout+len)) len++;
694 pfi->mnam = (char *)malloc(len+3);
695 if (pfi->mnam==NULL) goto err8;
696 strcpy (pfi->mnam,pout);
697
698 if (mflag) {
699 mfile = fopen (pout, "rb");
700 if (mfile==NULL) {
701 gaprnt (0,"Open Error: Can't open Station/Index map file ");
702 gaprnt (0,pout);
703 gaprnt (0,"\n");
704 goto err9;
705 }
706
707 if (mflag==2) goto skipread;
708
709 mflflg = 1;
710 swpflg = 0;
711
712 /* GRIB data - load the map INDEX data */
713
714 if (pfi->idxflg && pfi->type != 2) {
715
716 pindx = (struct gaindx *)malloc(sizeof(struct gaindx));
717 if (pindx==NULL) goto err8;
718 pfi->pindx = pindx;
719
720 /* check the version number */
721 fseek(mfile,1,0);
722 fread(&vermap,sizeof(unsigned char),1,mfile);
723 if(diag) printf("ddd gribmap vermap =%d\n",vermap);
724
725 if(vermap == 2) {
726
727 fseek(mfile,2,0);
728 fread(mrec,sizeof(unsigned char),4,mfile);
729 pindx->hinum=gagby(mrec,0,4);
730 if(diag) printf("ddd hinum = %d\n",pindx->hinum);
731
732 fread(mrec,sizeof(unsigned char),4,mfile);
733 pindx->hfnum=gagby(mrec,0,4);
734 if(diag) printf("ddd hfnum = %d\n",pindx->hfnum);
735
736 fread(mrec,sizeof(unsigned char),4,mfile);
737 pindx->intnum=gagby(mrec,0,4);
738 if(diag) printf("ddd intnum = %d\n",pindx->intnum);
739
740 fread(mrec,sizeof(unsigned char),4,mfile);
741 pindx->fltnum=gagby(mrec,0,4);
742 if(diag) printf("ddd fltnum = %d\n",pindx->fltnum);
743
744 /* skip the begining time struct info */
745 fread(mrec,sizeof(unsigned char),7,mfile);
746
747 pindx->hipnt = NULL;
748 pindx->hfpnt = NULL;
749 pindx->intpnt = NULL;
750 pindx->fltpnt = NULL;
751
752 if (pindx->hinum>0) {
753 pindx->hipnt = (int *)malloc(sizeof(int)*pindx->hinum);
754 if (pindx->hipnt==NULL) goto err8;
755
756 for(i=0;i<pindx->hinum;i++) {
757 fread(mrec,sizeof(unsigned char),4,mfile);
758 idum=gagby(mrec,0,4);
759 if(gagbb(mrec,0,1)) idum=-idum;
760 *(pindx->hipnt+i)=idum;
761 if(diag) printf("ddd 1 i = %d pindx->hipnt = %d\n",i,idum);
762 }
763
764 }
765
766
767 if (pindx->hfnum>0) {
768 pindx->hfpnt = (float *)malloc(sizeof(float)*pindx->hfnum);
769 if (pindx->hfpnt==NULL) goto err8;
770 fread (pindx->hfpnt,sizeof(float),pindx->hfnum,mfile);
771 }
772
773 if (pindx->intnum>0) {
774 pindx->intpnt = (int *)malloc(sizeof(int)*pindx->intnum);
775 if (pindx->intpnt==NULL) goto err8;
776 for(i=0;i<pindx->intnum;i++) {
777 fread(mrec,sizeof(unsigned char),4,mfile);
778 idum=gagby(mrec,0,4);
779 if(gagbb(mrec,0,1)) idum=-gagbb(mrec,1,31);
780 *(pindx->intpnt+i)=idum;
781 if(diag) printf("ddd 2 i = %d pindx->intpnt = %d\n",i,idum);
782 }
783
784 }
785
786 if (pindx->fltnum>0) {
787 pindx->fltpnt = (float *)malloc(sizeof(float)*pindx->fltnum);
788 if (pindx->fltpnt==NULL) goto err8;
789
790 for(i=0;i<pindx->fltnum;i++) {
791 fread(urec,sizeof(unsigned char),4,mfile);
792 fdum=ibm2flt(urec);
793 *(pindx->fltpnt+i)=fdum;
794 if(diag) printf("ddd 3 i = %d pindx->fltpnt = %g\n",i,fdum);
795 }
796
797 }
798
799 } else {
800
801 fseek (mfile,0L,0);
802 fread (pindx,sizeof(struct gaindx),1,mfile);
803 if (pindx->type>>24 > 0) swpflg=1;
804 if (swpflg) gabswp((float *)pindx,5);
805 pindx->hipnt = NULL;
806 pindx->hfpnt = NULL;
807 pindx->intpnt = NULL;
808 pindx->fltpnt = NULL;
809 if (pindx->hinum>0) {
810 pindx->hipnt = (int *)malloc(sizeof(int)*pindx->hinum);
811 if (pindx->hipnt==NULL) goto err8;
812 fread (pindx->hipnt,sizeof(int),pindx->hinum,mfile);
813 if (swpflg) gabswp((float *)(pindx->hipnt),pindx->hinum);
814 }
815 if (pindx->hfnum>0) {
816 pindx->hfpnt = (float *)malloc(sizeof(float)*pindx->hfnum);
817 if (pindx->hfpnt==NULL) goto err8;
818 fread (pindx->hfpnt,sizeof(float),pindx->hfnum,mfile);
819 if (swpflg) gabswp(pindx->hfpnt,pindx->hfnum);
820 }
821 if (pindx->intnum>0) {
822 pindx->intpnt = (int *)malloc(sizeof(int)*pindx->intnum);
823 if (pindx->intpnt==NULL) goto err8;
824 fread (pindx->intpnt,sizeof(int),pindx->intnum,mfile);
825 if (swpflg) gabswp((float *)(pindx->intpnt),pindx->intnum);
826 }
827 if (pindx->fltnum>0) {
828 pindx->fltpnt = (float *)malloc(sizeof(float)*pindx->fltnum);
829 if (pindx->fltpnt==NULL) goto err8;
830 fread (pindx->fltpnt,sizeof(float),pindx->fltnum,mfile);
831 if (swpflg) gabswp(pindx->fltpnt,pindx->fltnum);
832 }
833
834 }
835
836 /* map file processing -- stnmap - current NONINDEXED map (idxflg = 0) */
837
838 } else {
839
840 /*mf version detection mf*/
841
842 fread(rec,1,16,mfile); /* minimum map file is 16 bytes */
843
844 vermap=1;
845 if(strncmp(rec,"GrADS_stnmapV002",16) == 0) {
846 vermap=2;
847 }
848
849 if(vermap == 2) {
850
851 fread(mrec,sizeof(unsigned char),4,mfile);
852 idum=gagby(mrec,0,4);
853 if(gagbb(mrec,0,1)) idum=-idum;
854 mcnt=idum;
855
856 fread(mrec,sizeof(unsigned char),4,mfile);
857 idum=gagby(mrec,0,4);
858 if(gagbb(mrec,0,1)) idum=-idum;
859 maxlv=idum;
860
861 pfi->tstrt = (int *)malloc(sizeof(int)*mcnt);
862 pfi->tcnt = (int *)malloc(sizeof(int)*mcnt);
863
864 for(i=0;i<mcnt;i++) {
865 fread(mrec,sizeof(unsigned char),4,mfile);
866 idum=gagby(mrec,0,4);
867 if(gagbb(mrec,0,1)) idum=-idum;
868 *(pfi->tstrt+i)=idum;
869 }
870
871 for(i=0;i<mcnt;i++) {
872 fread(mrec,sizeof(unsigned char),4,mfile);
873 idum=gagby(mrec,0,4);
874 if(gagbb(mrec,0,1)) idum=-idum;
875 *(pfi->tcnt+i)=idum;
876 }
877
878
879 } else if(vermap ==1) {
880
881 /* reposition and read two local ints for version 1 map*/
882
883 fseek (mfile,0L,0);
884
885 fread (&mcnt,sizeof(int),1,mfile);
886 fread (&maxlv,sizeof(int),1,mfile);
887
888 if (maxlv>>24 > 0) swpflg = 1;
889 if (swpflg) {
890 gabswp((float *)(&mcnt),1);
891 gabswp((float *)(&maxlv),1);
892 }
893 pfi->tstrt = (int *)malloc(sizeof(int)*mcnt);
894 pfi->tcnt = (int *)malloc(sizeof(int)*mcnt);
895 fread (pfi->tstrt,sizeof(int),mcnt,mfile);
896 fread (pfi->tcnt,sizeof(int),mcnt,mfile);
897 if (swpflg) {
898 gabswp((float *)pfi->tstrt,1);
899 gabswp((float *)pfi->tcnt,1);
900 }
901 pfi->mtype = 1;
902
903 }
904
905 }
906 skipread:
907 fclose (mfile);
908 mflflg = 0;
909
910 } /* END OF mpfile check */
911
912
913 } else if (cmpwrd("toff",rec)) {
914 ch = nxtwrd(rec);
915 if (ch==NULL) {
916 gaprnt (0,"Open Error: Missing toff value\n");
917 goto err9;
918 }
919 pos = intprs(ch,&(pfi->tlpst));
920 if (pos==NULL || pfi->tlpst>=pfi->dnum[3]) {
921 gaprnt (0,"Open Error: Invalid toff value\n");
922 goto err9;
923 }
924 pfi->tlpflg = 1;
925
926
927 } else if (cmpwrd("undef",rec)) {
928 ch = nxtwrd(mrec);
929 if (ch==NULL) {
930 gaprnt (0,"Open Error: Missing undef value\n");
931 goto err9;
932 }
933
934 pos = valprs(ch,&(pfi->undef));
935 if (pos==NULL) {
936 gaprnt (0,"Open Error: Invalid undef value\n");
937 goto err9;
938 }
939
940 /* Get the undef attribute name, if it's there */
941 if ( (ch=nxtwrd(ch))!=NULL ) {
942 len = 0;
943 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
944 pfi->undefattr = (char *)malloc(len+3);
945 for (i=0; i<len; i++) *(pfi->undefattr+i) = *(ch+i);
946 *(pfi->undefattr+len) = '\0';
947 /* Set the undef attribute flag */
948 pfi->undefattrflg = 1;
949 }
950
951 if (SETMISS) {
952 pfi->ulow = fabs(pfi->undef/EPSILON);
953 pfi->uhi = pfi->undef + pfi->ulow;
954 pfi->ulow = pfi->undef - pfi->ulow;
955 }
956 flgs[4] = 0;
957
958 } else if (cmpwrd("pdef",rec)) {
959 if ( (ch = nxtwrd(rec)) == NULL) goto errm;
960 /* parse the i and j dimensions of the pre-projected grid */
961 if ( (pos = intprs(ch,&(pfi->ppisiz)))==NULL) goto errm;
962 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
963 if ( (pos = intprs(ch,&(pfi->ppjsiz)))==NULL) goto errm;
964 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
965 /* set the pre-projected grid type and wind rotation flags */
966 if (cmpwrd("nps",ch)) {pfi->ppflag=1; pfi->ppwrot=1; cnt=4;}
967 else if (cmpwrd("sps",ch)) {pfi->ppflag=2; pfi->ppwrot=1; cnt=4;}
968 else if (cmpwrd("lcc",ch)) {pfi->ppflag=3; pfi->ppwrot=0; cnt=9;}
969 else if (cmpwrd("lccr",ch)) {pfi->ppflag=3; pfi->ppwrot=1; cnt=9;}
970 else if (cmpwrd("eta.u",ch)) {pfi->ppflag=4; pfi->ppwrot=1; cnt=4;}
971 else if (cmpwrd("pse",ch)) {pfi->ppflag=5; pfi->ppwrot=0; cnt=7;}
972 else if (cmpwrd("ops",ch)) {pfi->ppflag=6; pfi->ppwrot=0; cnt=8;}
973 else if (cmpwrd("bilin",ch)) {pfi->ppflag=7; pfi->ppwrot=1; cnt=0;}
974 else if (cmpwrd("file",ch)) {pfi->ppflag=8; pfi->ppwrot=1; cnt=1;}
975 else goto errm;
976 /* parse the pre-projected grid parameters */
977 for (i=0; i<cnt; i++) {
978 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
979 if ( (pos = valprs(ch,&(pfi->ppvals[i])))==NULL) goto errm;
980 }
981 /* check "num" argument to pdef file option */
982 if (pfi->ppflag==8) {
983 i = (int)(pfi->ppvals[0]+0.1);
984 if (i<1 || i>9) goto errm;
985 }
986 /* parse file type, byte order, and name for pdef 'bilin' and 'file' */
987 if (pfi->ppflag==7 || pfi->ppflag==8) {
988 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
989 if (cmpwrd("stream",ch)) pdefop1 = 1;
990 else if (cmpwrd("sequential",ch)) pdefop1 = 2;
991 else goto errm;
992
993 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
994 if (cmpwrd("binary",ch)) pdefop2 = 1;
995 else if (cmpwrd("binary-big",ch)) pdefop2 = 2;
996 else if (cmpwrd("binary-little",ch)) pdefop2 = 3;
997 else if (cmpwrd("packed",ch)) pdefop2 = 4;
998 else goto errm;
999
1000 if ( (ch = nxtwrd(ch)) == NULL) goto errm;
1001 ch = mrec + (ch-rec);
1002 if (*ch=='^' || *ch=='$') {
1003 fnmexp (pdefnm,ch,name);
1004 } else {
1005 getwrd (pdefnm,ch,256);
1006 }
1007 /* open the pdef file */
1008 pdfi = fopen(pdefnm,"rb");
1009 if (pdfi==NULL) {
1010 sprintf (pout, " Error opening pdef file: %s\n",pdefnm);
1011 gaprnt (0,pout);
1012 goto errm;
1013 }
1014 }
1015
1016 } else if (cmpwrd("vectorpairs",rec)) {
1017 if ( (ch=nxtwrd(mrec))==NULL ) {
1018 gaprnt (1,"Description file warning: No vector pairs listed\n");
1019 } else {
1020 if ((vectorpairs = (char *)malloc(strlen(ch)+1)) == NULL) {
1021 gaprnt(0,"malloc error for list of vector pairs\n");
1022 return(1);
1023 }
1024 getstr(vectorpairs,ch,strlen(ch)+1);
1025 }
1026
1027 } else if (cmpwrd("xdef",rec)) {
1028
1029 if (pfi->type == 2) continue;
1030 if ( (ch = nxtwrd(rec)) == NULL) goto err1;
1031 if ( (pos = intprs(ch,&(pfi->dnum[0])))==NULL) goto err1;
1032 if (*pos!=' ') goto err1;
1033 if ( (ch = nxtwrd(ch))==NULL) goto err2;
1034 if (cmpwrd("linear",ch)) {
1035 rc = deflin(ch, pfi, 0, 0);
1036 if (rc==-1) goto err8;
1037 if (rc) goto err9;
1038 v2 = *(pfi->grvals[0]);
1039 v1 = *(pfi->grvals[0]+1) + v2;
1040 temp = v1+((float)(pfi->dnum[0]))*v2;
1041 temp=temp-360.0;
1042 if (fabs(temp-v1)<0.01) pfi->wrap = 1;
1043 }
1044 else if (cmpwrd("levels",ch)) {
1045 if (pfi->dnum[0]<1) goto err7;
1046 rc = deflev (ch, rec, pfi, 0);
1047 if (rc==-1) goto err8;
1048 if (rc) goto err9;
1049 } else goto err2;
1050 flgs[0] = 0;
1051
1052 } else if (cmpwrd("ydef",rec)) {
1053 if (pfi->type == 2) continue;
1054 if ( (ch = nxtwrd(rec)) == NULL) goto err1;
1055 if ( (pos = intprs(ch,&(pfi->dnum[1])))==NULL) goto err1;
1056 if (*pos!=' ') goto err1;
1057 if ( (ch = nxtwrd(ch))==NULL) goto err2;
1058 if (cmpwrd("linear",ch)) {
1059 rc = deflin(ch, pfi, 1, 0);
1060 if (rc==-1) goto err8;
1061 if (rc) goto err9;
1062 } else if (cmpwrd("levels",ch)) {
1063 if (pfi->dnum[1]<1) goto err7;
1064 rc = deflev (ch, rec, pfi, 1);
1065 if (rc==-1) goto err8;
1066 if (rc) goto err9;
1067 } else if (cmpwrd("gausr40",ch)) {
1068 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1069 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1070 pfi->grvals[1] = gagaus(i,pfi->dnum[1]);
1071 if (pfi->grvals[1]==NULL) goto err9;
1072 pfi->abvals[1] = pfi->grvals[1];
1073 pfi->ab2gr[1] = lev2gr;
1074 pfi->gr2ab[1] = gr2lev;
1075 pfi->linear[1] = 0;
1076 } else if (cmpwrd("mom32",ch)) {
1077 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1078 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1079 pfi->grvals[1] = gamo32(i,pfi->dnum[1]);
1080 if (pfi->grvals[1]==NULL) goto err9;
1081 pfi->abvals[1] = pfi->grvals[1];
1082 pfi->ab2gr[1] = lev2gr;
1083 pfi->gr2ab[1] = gr2lev;
1084 pfi->linear[1] = 0;
1085 } else if (cmpwrd("gaust62",ch)) {
1086 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1087 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1088 pfi->grvals[1] = gagst62(i,pfi->dnum[1]);
1089 if (pfi->grvals[1]==NULL) goto err9;
1090 pfi->abvals[1] = pfi->grvals[1];
1091 pfi->ab2gr[1] = lev2gr;
1092 pfi->gr2ab[1] = gr2lev;
1093 pfi->linear[1] = 0;
1094 } else if (cmpwrd("gausr30",ch)) {
1095 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1096 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1097 pfi->grvals[1] = gags30(i,pfi->dnum[1]);
1098 if (pfi->grvals[1]==NULL) goto err9;
1099 pfi->abvals[1] = pfi->grvals[1];
1100 pfi->ab2gr[1] = lev2gr;
1101 pfi->gr2ab[1] = gr2lev;
1102 pfi->linear[1] = 0;
1103 } else if (cmpwrd("gausr20",ch)) {
1104 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1105 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1106 pfi->grvals[1] = gags20(i,pfi->dnum[1]);
1107 if (pfi->grvals[1]==NULL) goto err9;
1108 pfi->abvals[1] = pfi->grvals[1];
1109 pfi->ab2gr[1] = lev2gr;
1110 pfi->gr2ab[1] = gr2lev;
1111 pfi->linear[1] = 0;
1112 } else if (cmpwrd("gausr15",ch)) {
1113 if ( (ch = nxtwrd(ch))==NULL) goto err3;
1114 if ( (pos = intprs(ch,&i))==NULL) goto err3;
1115 pfi->grvals[1] = gags15(i,pfi->dnum[1]);
1116 if (pfi->grvals[1]==NULL) goto err9;
1117 pfi->abvals[1] = pfi->grvals[1];
1118 pfi->ab2gr[1] = lev2gr;
1119 pfi->gr2ab[1] = gr2lev;
1120 pfi->linear[1] = 0;
1121 } else goto err2;
1122 flgs[1] = 0;
1123
1124 } else if (cmpwrd("zdef",rec)) {
1125 if (pfi->type == 2) continue;
1126 if ( (ch = nxtwrd(rec)) == NULL) goto err1;
1127 if ( (pos = intprs(ch,&(pfi->dnum[2])))==NULL) goto err1;
1128 if (*pos!=' ') goto err1;
1129 if ( (ch = nxtwrd(ch))==NULL) goto err2;
1130 if (cmpwrd("linear",ch)) {
1131 rc = deflin(ch, pfi, 2, 0);
1132 if (rc==-1) goto err8;
1133 if (rc) goto err9;
1134 }
1135 else if (cmpwrd("levels",ch)) {
1136 if (pfi->dnum[2]<1) goto err7;
1137 rc = deflev (ch, rec, pfi, 2);
1138 if (rc==-1) goto err8;
1139 if (rc) goto err9;
1140 } else goto err2;
1141 flgs[2] = 0;
1142
1143 } else if (cmpwrd("tdef",rec)) {
1144 if ( (ch = nxtwrd(rec)) == NULL) goto err1;
1145 if ( (pos = intprs(ch,&(pfi->dnum[3])))==NULL) goto err1;
1146 if (*pos!=' ') goto err1;
1147 if ( (ch = nxtwrd(ch))==NULL) goto err2;
1148 if (cmpwrd("linear",ch)) {
1149 if ( (ch = nxtwrd(ch))==NULL) goto err3a_tdef;
1150 tdef.yr = -1000;
1151 tdef.mo = -1000;
1152 tdef.dy = -1000;
1153 if ( (pos = adtprs(ch,&tdef,&dt1))==NULL) goto err3b_tdef;
1154 if (*pos!=' ' || dt1.yr == -1000 || dt1.mo == -1000.0 ||
1155 dt1.dy == -1000) goto err3c_tdef;
1156 if ( (ch = nxtwrd(ch))==NULL) goto err4a_tdef;
1157 if ( (pos = rdtprs(ch,&dt2))==NULL) goto err4b_tdef;
1158 v1 = (dt2.yr * 12) + dt2.mo;
1159 v2 = (dt2.dy * 1440) + (dt2.hr * 60) + dt2.mn;
1160 /*mf --- check if 0 dt ---mf*/
1161 if( (v1 == 0) && (v2 == 0) ) goto err4c_tdef;
1162 vals = (float *)malloc(sizeof(float)*8);
1163 if (vals==NULL) goto err8;
1164 *(vals) = dt1.yr;
1165 *(vals+1) = dt1.mo;
1166 *(vals+2) = dt1.dy;
1167 *(vals+3) = dt1.hr;
1168 *(vals+4) = dt1.mn;
1169 *(vals+5) = v1;
1170 *(vals+6) = v2;
1171 *(vals+7) = -999.9;
1172 pfi->grvals[3] = vals;
1173 pfi->abvals[3] = vals;
1174 pfi->linear[3] = 1;
1175 } else goto err2;
1176 flgs[3] = 0;
1177
1178 } else if (cmpwrd("vars",rec)) {
1179 if ( (ch = nxtwrd(rec)) == NULL) goto err5;
1180 if ( (pos = intprs(ch,&(pfi->vnum)))==NULL) goto err5;
1181 size = pfi->vnum * (sizeof(struct gavar) + 7 );
1182 pvar = (struct gavar *)malloc(size);
1183 pfi->pvar1 = pvar;
1184 i = 0;
1185 while (i<pfi->vnum) {
1186 if (fgets(rec,512,descr)==NULL) {
1187 gaprnt (0,"Open Error: Unexpected EOF reading variables\n");
1188 sprintf (pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
1189 gaprnt (2,pout);
1190 goto retrn;
1191 }
1192
1193 /* Remove any leading blanks from rec */
1194 reclen = strlen(rec);
1195 jj = 0;
1196 while (jj<reclen && rec[0]==' ') {
1197 for (ii=0; ii<reclen; ii++) rec[ii] = rec[ii+1];
1198 jj++;
1199 }
1200
1201 /* Keep mixed case and lower case versions of rec handy */
1202 strcpy (mrec,rec);
1203 lowcas(rec);
1204
1205 /* Allow comments between VARS and ENDVARS */
1206 if ( (strncmp(mrec,"*",1)==0) || (strncmp(mrec,"#",1)==0) || !isalnum(*(mrec)) ) {
1207 /* Parse comment if it contains attribute metadata */
1208 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0)) {
1209 rc = ddfattr (mrec, pfi);
1210 if (rc == -1) goto retrn;
1211 else continue;
1212 }
1213 else continue;
1214 }
1215 if (cmpwrd("endvars",rec)) {
1216 gaprnt (0,"Open Error: Unexpected ENDVARS record\n");
1217 sprintf (pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
1218 gaprnt (2,pout);
1219 goto err9;
1220 }
1221
1222 /* get abbrv and full variable name if there */
1223 pvar->longnm = NULL;
1224 if ((getvnm(pvar, mrec))!=0) goto err6;
1225
1226 /* get the number of vertical levels */
1227 if ( (ch=nxtwrd(rec))==NULL) goto err6;
1228 if ( (pos=intprs(ch,&(pvar->levels)))==NULL ) goto err6;
1229 if ( (ch=nxtwrd(ch))==NULL) goto err6;
1230
1231 /* parse the units field */
1232 for (j=0;j<16;j++) pvar->units[j] = -999;
1233 j = 0;
1234 while (1) {
1235 if (*ch=='x'||*ch=='y'||*ch=='z'||*ch=='t') {
1236 if (*(ch+1)!=',' && *(ch+1)!=' ') goto err6;
1237 if (*ch=='x') pvar->units[j] = -100;
1238 if (*ch=='y') pvar->units[j] = -101;
1239 if (*ch=='z') pvar->units[j] = -102;
1240 if (*ch=='t') pvar->units[j] = -103;
1241 ch++;
1242 } else {
1243 if ( (ch=intprs(ch,&(pvar->units[j])))==NULL ) goto err6;
1244 /* no negative array indices for ncflag files */
1245 if ((pfi->ncflg) && (pvar->units[j] < 0)) goto err6;
1246 }
1247 while (*ch==' ') ch++;
1248 if (*ch=='\0' || *ch=='\n') goto err6;
1249 if (*ch!=',') break;
1250 ch++;
1251 while (*ch==' ') ch++;
1252 if (*ch=='\0' || *ch=='\n') goto err6;
1253 j++;
1254 if (j>15) goto err6;
1255 }
1256
1257 /* parse the variable description */
1258 getstr (pvar->varnm,mrec+(ch-rec),127);
1259
1260 /* initialize the vector pair number and the "is u" flag*/
1261 pvar->vecpair = -999;
1262 pvar->isu = 0;
1263
1264 /* initialize the var_z counter for NASA GLA format */
1265 pvar->var_z = 1;
1266 if(pvar->units[0] == -1 && pvar->units[1] == 10 ) {
1267 nzstride++;
1268 }
1269
1270 /* var_t is for DRS var-t transforms */
1271 pvar->var_t = 0;
1272 if(pvar->units[0] == -1 && pvar->units[1] == 20 ) pvar->var_t = 1;
1273
1274 /* x-y transpose for lat/lon vice lon/lat data VERY INEFFICIENT!!! */
1275 pvar->y_x = 0;
1276 if(pvar->units[0] == -1 && pvar->units[1] == 30) pvar->y_x = 1;
1277
1278 /*mf 961126 - integer data types mf*/
1279 pvar->dfrm = 0;
1280 if( (pvar->units[0] == -1 && pvar->units[1] == 40 ) &&
1281 (pvar->units[2] >=1 && pvar->units[2] <= 4 ) ) {
1282 pvar->dfrm= pvar->units[2];
1283 if ( pvar->units[2] == 2 ) pvar->dfrm=2;
1284 if ( pvar->units[2] == 2 && pvar->units[3] == -1) pvar->dfrm=-2;
1285 }
1286
1287 #if GRADS_CRAY == 1
1288 /* 32-bit big endian ieee to cray float */
1289 if( (pvar->units[0]==-1 && pvar->units[1]==50) || crayflg ) pvar->dfrm=8;
1290 #endif
1291 i++; pvar++;
1292 }
1293
1294 /* Get ENDVARS statement and any additional comments */
1295 if (fgets(rec,512,descr)==NULL) {
1296 gaprnt (0,"Open Error: Missing ENDVARS statement.\n");
1297 goto retrn;
1298 }
1299 /* Remove any leading blanks from rec */
1300 reclen = strlen(rec);
1301 jj = 0;
1302 while (jj<reclen && rec[0]==' ') {
1303 for (ii=0; ii<reclen; ii++) rec[ii] = rec[ii+1];
1304 jj++;
1305 }
1306 /* Keep mixed case and lower case versions handy */
1307 strcpy (mrec,rec);
1308 lowcas(rec);
1309 while (!cmpwrd("endvars",rec)) {
1310 /* see if it's an attribute comment, otherwise send error message */
1311 if ((strncmp("*:attr",mrec,6)==0) || (strncmp("@",mrec,1)==0)) {
1312 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
1313 }
1314 else {
1315 sprintf(pout,"Open Error: Looking for \"endvars\", found \"%s\" instead.\n",rec);
1316 gaprnt (0,pout);
1317 goto err9;
1318 }
1319 if (fgets(rec,512,descr)==NULL) {
1320 gaprnt (0,"Open Error: Missing ENDVARS statement.\n");
1321 goto retrn;
1322 }
1323 }
1324
1325 flgs[6] = 0;
1326 } else {
1327 /* parse error of .ctl file */
1328 gaprnt (0,"Open Error: Unknown keyword in description file\n");
1329 goto err9;
1330 }
1331 }
1332
1333 err=0;
1334 for (i=0; i<7; i++) {
1335 if (flgs[i]) {
1336 sprintf (pout,"Open Error: missing %s record \n",errs[i]);
1337 gaprnt (0,pout);
1338 err=1;
1339 }
1340 }
1341
1342 if (err) goto retrn;
1343
1344
1345 /* Done scanning. Check if scanned stuff makes sense,
1346 and then set things up correctly */
1347
1348 npairs = 0;
1349 if (vectorpairs) {
1350 /* parse the vector pairs */
1351 if ((pair = (char *)malloc(strlen(vectorpairs)+1)) == NULL) {
1352 gaprnt(0,"malloc error for list of vector pairs\n");
1353 err=1;
1354 }
1355 else {
1356 while (1) {
1357 /* copy the comma-delimited pair of variable names */
1358 getwrd(pair,vectorpairs,strlen(vectorpairs));
1359 i=0;
1360 while (1) {
1361 if (*(pair+i)==',') {
1362 /* get the two variable names that comprise the vector pair */
1363 getstr(var1, pair, i);
1364 getstr(var2, pair+(i+1), strlen(pair)-i+1);
1365 npairs++;
1366
1367 /* loop through variable list */
1368 foundvar1=foundvar2=0;
1369 pvar = pfi->pvar1;
1370 for (j=1; j<=pfi->vnum; j++) {
1371 /* if variable names match, set flags */
1372 if (cmpwrd(pvar->abbrv,var1)) foundvar1=j;
1373 if (cmpwrd(pvar->abbrv,var2)) foundvar2=j;
1374 pvar++;
1375 }
1376 if (foundvar1 && foundvar2) {
1377 pvar = pfi->pvar1;
1378 for (j=1; j<=pfi->vnum; j++) {
1379 /* if we've found both variables, set pvar->vecpair */
1380 if (cmpwrd(pvar->abbrv,var1)) {
1381 pvar->vecpair=npairs;
1382 pvar->isu=1; /* trip flag for u-component */
1383 }
1384 if (cmpwrd(pvar->abbrv,var2)) {
1385 pvar->vecpair=npairs;
1386 }
1387 pvar++;
1388 }
1389 }
1390 else {
1391 sprintf(pout,"Warning: VECTORPAIRS variables %s,%s were not found \n",var1,var2);
1392 gaprnt(1,pout);
1393 }
1394 break;
1395 }
1396 i++;
1397 }
1398 /* move pointer forward one word */
1399 if ((vectorpairs = nxtwrd (vectorpairs)) == NULL) break;
1400 }
1401 }
1402 free(pair);
1403 free(vectorpairs);
1404 }
1405
1406 /* Find variables with units[0]=33,34 -- vector pairs that havn't already been flagged */
1407 pvar=pfi->pvar1;
1408 for (j=1; j<=pfi->vnum; j++) {
1409 /* Look for a variable with units[0]==33,34 that hasn't been handled yet */
1410 if ((pvar->units[0]==33 || pvar->units[0]==34) && (pvar->vecpair<0)) {
1411 if (pvar->units[0]==33) { rc = 34; }
1412 else { rc = 33; }
1413 /* Look for a matching variable with all units fields equal */
1414 pvar2 = (struct gavar *)malloc(size);
1415 if (pvar2==NULL) {
1416 gaprnt (0,"Memory allocation error for pvar2 in gaddes.c \n");
1417 err=1;
1418 goto retrn;
1419 }
1420 pvar2 = pfi->pvar1;
1421 i = 0;
1422 while (i<pfi->vnum) {
1423 if (pvar2->vecpair<0 &&
1424 pvar2->units[0]==rc &&
1425 (pvar->units[1]==pvar2->units[1]) &&
1426 (pvar->units[2]==pvar2->units[2]) &&
1427 (pvar->units[3]==pvar2->units[3]) ) break;
1428 pvar2++; i++;
1429 }
1430 if (i<pfi->vnum) { /* We've got a match! */
1431 npairs++;
1432 /* set the gavar parameters */
1433 pvar->vecpair=npairs;
1434 pvar2->vecpair=npairs;
1435 if (pvar->units[0]==33) { pvar->isu=1; }
1436 else { pvar2->isu=1; }
1437 }
1438 }
1439 pvar++;
1440 }
1441 if (err) goto retrn; /* end of vector pair management */
1442
1443 if (pfi->type>1 && mflag==1) {
1444 if (mcnt==-1) {
1445 gaprnt (0,"Open Error: missing STNMAP record\n");
1446 err=1;
1447 } else if (mcnt != pfi->dnum[3]) {
1448 gaprnt (0,"Open Error: Inconsistent time count\n");
1449 sprintf (pout," Count in station map file = %i\n",mcnt);
1450 gaprnt (0,pout);
1451 sprintf (pout," Count in descriptor file = %i\n",pfi->dnum[3]);
1452 gaprnt (0,pout);
1453 err=1;
1454 }
1455 }
1456
1457 if (err) goto retrn;
1458
1459 /* Figure out locations of variables within a time group */
1460 pvar = pfi->pvar1;
1461
1462 /* Grid data */
1463 if (pfi->type==1) {
1464 pfi->gsiz = pfi->dnum[0] * pfi->dnum[1];
1465 if (pfi->ppflag) pfi->gsiz = pfi->ppisiz * pfi->ppjsiz;
1466
1467 /* add a constant xy header */
1468 if (pfi->xyhdr) {
1469 pfi->gsiz = pfi->gsiz + pfi->xyhdr;
1470 }
1471
1472 if (pfi->seqflg) {
1473 pfi->gsiz+=2;
1474 if (hdrb>0) hdrb+=2;
1475 pvar->offset = 1+hdrb;
1476 acum = 1+hdrb;
1477 } else {
1478 pvar->offset = hdrb;
1479 acum = hdrb;
1480 }
1481
1482 levs = pvar->levels;
1483 if (levs==0) levs=1;
1484 pvar->recoff = 0;
1485 recacm = 0;
1486 pvar++;
1487
1488 acumvz=acum;
1489
1490 for (i=1; i<pfi->vnum; i++) {
1491
1492 /* NASA GLA FORMAT CHECKS */
1493 /* upper air fields which var and z are transposed*/
1494 if((pvar->units[0]==-1) &&
1495 (pvar->units[1]==10) &&
1496 (pvar->units[2]==1) ) {
1497
1498 levsua = pvar->levels;
1499 acum = acum + pfi->gsiz;
1500 pvar->var_z = nzstride;
1501
1502 /* diagnstotic fields AFTER upper air fields which are in GrADS normal order */
1503
1504 } else if((pvar->units[0]==-1) &&
1505 (pvar->units[1]==10) &&
1506 (pvar->units[2]==2) ) {
1507 if(first) {
1508 acumstride = acumstride + nzstride*levsua*pfi->gsiz;
1509 acum = acumstride + (pvar->levels*pfi->gsiz) ;
1510 first = 0 ;
1511 } else {
1512 acum = acum + (levs*pfi->gsiz);
1513 }
1514
1515 } else if(pvar->var_t) { /* DRS transposition of var and t */
1516
1517 /* time template, read the third unit param for the size of each chunk in a file */
1518 if(pfi->tmplat) {
1519 if(pvar->units[2] != -999) {
1520 acum = acum + levs*(pfi->gsiz)*(pvar->units[1]);
1521 } else {
1522 gaprnt (0,"Using time templat and 4-D variables, # times / file not specified\n");
1523 gaprnt (0,"Defaulting to the number of times in the .ctl file\n");
1524 acum = acum + levs*(pfi->gsiz)*(pfi->dnum[3]);
1525 }
1526 } else {
1527 acum = acum + levs*(pfi->gsiz)*(pfi->dnum[3]);
1528 }
1529 } else { /* simple GrADS */
1530 acum = acum + (levs*pfi->gsiz);
1531 acumstride = acum ;
1532 }
1533 recacm += levs;
1534 pvar->offset = acum;
1535 pvar->recoff = recacm;
1536 levs = pvar->levels;
1537 if (levs==0) levs=1;
1538 pvar++;
1539 }
1540
1541 recacm += levs;
1542
1543 /*mf 960514 correct for case where the last variable is a NASA UA */
1544 pvar--;
1545 if( (pvar->units[0]==-1) &&
1546 (pvar->units[1]==10) &&
1547 (pvar->units[2]==1) ) {
1548 acum = acumvz + recacm*pfi->gsiz;
1549 }
1550 else {
1551 acum = acum + (levs*pfi->gsiz);
1552 }
1553
1554 /* time chunk header; the default is 0 */
1555 if(pfi->seqflg && pfi->thdr>0) {
1556 pfi->thdr+=2;
1557 }
1558 pfi->tsiz = acum + pfi->thdr;
1559 pfi->trecs = recacm;
1560 if (pfi->seqflg) pfi->tsiz-=1;
1561 pfi->tsiz += trlb;
1562
1563 } else { /* non grid data */
1564
1565 for (i=0; i<pfi->vnum; i++) {
1566 if (pvar->levels!=0) break;
1567 pvar->offset = i;
1568 pvar++;
1569 }
1570 for (j=i; j<pfi->vnum; j++) {
1571 if (pvar->levels==0) {
1572 if (!(pfi->bufrflg)) { /* order of variables doesn't matter for BUFR data */
1573 gaprnt (0,"Open Error: Variables out of order\n");
1574 gaprnt (0," Non-vertical variables must go first\n");
1575 goto retrn;
1576 }
1577 }
1578 pvar->offset = j-i;
1579 pvar++;
1580 }
1581 pfi->lvnum = pfi->vnum - i;
1582 pfi->ivnum = i;
1583 }
1584
1585 /* set the global calendar and check if we are trying to change with a new file...
1586 we do this here to set the calandar for templating */
1587
1588 if(mfcmn.cal365<0) {
1589 mfcmn.cal365=pfi->calendar;
1590 } else {
1591 if(pfi->calendar != mfcmn.cal365) {
1592 gaprnt(0,"Attempt to change the global calendar...\n");
1593 if(mfcmn.cal365) {
1594 gaprnt(0,"The calendar is NOW 365 DAYS and you attempted to open a standard calendar file\n");
1595 } else {
1596 gaprnt(0,"The calendar is NOW STANDARD and you attempted to open a 365-day calendar file\n");
1597 }
1598 goto retrn;
1599 }
1600 }
1601
1602 /* Allocate an I/O buffer the size of one row */
1603 if (pfi->type > 1) {
1604 if (pfi->bufrflg) maxlv=1; /* Not used for BUFR interface, set to 1 */
1605 size = maxlv * sizeof(float);
1606 } else {
1607 size = pfi->dnum[0] * sizeof(float);
1608 }
1609 pfi->rbuf = (float *)malloc(size);
1610 if (pfi->idxflg) pfi->pbuf = (char *)malloc(size);
1611
1612 /* If a pre-projected grid, set up the interpolation
1613 constants. If we're just checking the descriptor (mflag==2),
1614 no need to do this */
1615
1616 if (pfi->ppflag && mflag!=2) {
1617 rc = gappcn(pfi,pdefop1,pdefop2);
1618 if (rc) goto retrn;
1619 }
1620
1621 /* If the file name is a time series template, figure out
1622 which times go with which files, so we don't waste a lot
1623 of time later opening and closing files unnecessarily. */
1624
1625 /* BUFR files are treated like templated files, so that the
1626 data file isn't parsed until an I/O request is made */
1627
1628 if (pfi->tmplat || pfi->bufrflg==1) {
1629 pfi->fnums = (int *)malloc(sizeof(int)*pfi->dnum[3]);
1630 if (pfi->fnums==NULL) goto err8;
1631 j = 1;
1632 gr2t(pfi->grvals[3],1.0,&tdefi);
1633 ch = gafndt(pfi->name,&tdefi,&tdefi,pfi->abvals[3],pfi->pchsub1,1);
1634 if (ch==NULL) goto err8;
1635 *(pfi->fnums) = j;
1636 for (i=2; i<=pfi->dnum[3]; i++) {
1637 gr2t(pfi->grvals[3],(float)i,&tdef);
1638 pos = gafndt(pfi->name,&tdef,&tdefi,pfi->abvals[3],pfi->pchsub1,i);
1639 if (pos==NULL) goto err8;
1640 if (strcmp(ch,pos)!=0) {
1641 j = i;
1642 free(ch);
1643 ch = pos;
1644 }
1645 *(pfi->fnums+i-1) = j;
1646 }
1647 free(ch);
1648 pfi->fnumc = 0;
1649 }
1650
1651 fclose (descr);
1652 if (pdfi) fclose(pdfi);
1653
1654 return(0);
1655
1656 errm:
1657 gaprnt(0,"Open Error: Invalid pdef record.\n");
1658 pfi->ppflag = 0;
1659 goto err9;
1660
1661 err1:
1662 gaprnt (0,"Open Error: Missing or invalid dimension size.\n");
1663 goto err9;
1664
1665 err2:
1666 gaprnt (0,"Open Error: Missing or invalid dimension");
1667 gaprnt (0," scaling type\n");
1668 goto err9;
1669
1670 err3a_tdef:
1671 gaprnt (0,"Open Error: Start Time missing in tdef card");
1672 gaprnt (0," starting value\n");
1673 goto err9;
1674
1675 err3b_tdef:
1676 gaprnt (0,"Open Error: Invalid start time in tdef card");
1677 gaprnt (0," starting value\n");
1678 goto err9;
1679
1680 err3c_tdef:
1681 gaprnt (0,"Open Error: Missing or invalid dimension");
1682 gaprnt (0," starting value\n");
1683 goto err9;
1684
1685 err3:
1686 gaprnt (0,"Open Error: Missing or invalid dimension");
1687 gaprnt (0," starting value\n");
1688 goto err9;
1689
1690 err4a_tdef:
1691 gaprnt (0,"Open Error: Time increment missing in tdef\n");
1692 gaprnt (0," use 1 for single time data\n");
1693 goto err9;
1694
1695 err4b_tdef:
1696 gaprnt (0,"Open Error: Invalid time increment in tdef\n");
1697 gaprnt (0," use 1 for single time data\n");
1698 goto err9;
1699
1700 err4c_tdef:
1701 gaprnt (0,"Open Error: 0 time increment in tdef\n");
1702 gaprnt (0," use 1 for single time data\n");
1703 goto err9;
1704
1705 err5:
1706 gaprnt (0,"Open Error: Missing or invalid variable");
1707 gaprnt (0," count\n");
1708 goto err9;
1709
1710 err6:
1711 gaprnt (0,"Open Error: Invalid variable record\n");
1712 goto err9;
1713
1714 err6a:
1715 gaprnt (0,"Open Error: Invalid x,y pair\n");
1716 goto err9;
1717
1718 err7:
1719 gaprnt (0,"Open Error: Invalid number of levels\n");
1720 goto err9;
1721
1722 err8:
1723 gaprnt (0,"Open Error: Memory allocation Error in gaddes.c\n");
1724 goto retrn;
1725
1726 err9:
1727 gaprnt (0," --> The invalid description file record is: \n");
1728 gaprnt (0," --> ");
1729 gaprnt (0,mrec);
1730 gaprnt (0,"\n");
1731
1732 retrn:
1733 gaprnt (0," The data file was not opened. \n");
1734 fclose (descr);
1735 if (mflflg) fclose(mfile);
1736 if (pdfi) fclose(pdfi);
1737 return(1);
1738
1739 }
1740
1741
1742
1743 /* Process linear scaling args */
1744
deflin(char * ch,struct gafile * pfi,int dim,int flag)1745 int deflin (char *ch, struct gafile *pfi, int dim, int flag) {
1746 float *vals,v1,v2;
1747
1748 vals = (float *)malloc(sizeof(float)*6);
1749 if (vals==NULL) return (-1);
1750
1751 if ( (ch = nxtwrd(ch))==NULL) goto err1;
1752 if ( valprs(ch,&v1)==NULL) goto err1;
1753 if (flag) v2 = 1.0;
1754 else {
1755 if ( (ch = nxtwrd(ch))==NULL) goto err2;
1756 if ( valprs(ch,&v2)==NULL) goto err2;
1757 }
1758 if (dim<3 && v2<=0.0) goto err2;
1759 *(vals+1) = v1 - v2;
1760 *(vals) = v2;
1761 *(vals+2) = -999.9;
1762 pfi->grvals[dim] = vals;
1763 *(vals+4) = -1.0 * ( (v1-v2)/v2 );
1764 *(vals+3) = 1.0/v2;
1765 *(vals+5) = -999.9;
1766 pfi->abvals[dim] = vals+3;
1767 pfi->ab2gr[dim] = liconv;
1768 pfi->gr2ab[dim] = liconv;
1769 pfi->linear[dim] = 1;
1770 return (0);
1771
1772 err1:
1773 gaprnt (0,"Open Error: Missing or invalid dimension");
1774 gaprnt (0," starting value\n");
1775 free (vals);
1776 return (1);
1777
1778 err2:
1779 gaprnt (0,"Open Error: Missing or invalid dimension");
1780 gaprnt (0," increment value\n");
1781 free (vals);
1782 return (1);
1783 }
1784
1785 /* Process levels values in def record */
1786 /* Return codes: -1 is memory allocation error, 1 is other error */
1787
deflev(char * ch,char * rec,struct gafile * pfi,int dim)1788 int deflev (char *ch, char *rec, struct gafile *pfi, int dim) {
1789 float *vvs,*vals,v1;
1790 int i;
1791
1792 if (pfi->dnum[dim]==1) {
1793 i = deflin (ch, pfi, dim, 1);
1794 return (i);
1795 }
1796
1797 vals = (float *)malloc((pfi->dnum[dim]+5)*sizeof(float));
1798 if (vals==NULL) return (-1);
1799
1800 vvs = vals;
1801 *vvs = (float)pfi->dnum[dim];
1802 vvs++;
1803 for (i=0; i<pfi->dnum[dim]; i++) {
1804 if ( (ch = nxtwrd(ch))==NULL) {
1805 if (fgets(rec,256,descr)==NULL) goto err2;
1806 ch = rec;
1807 while (*ch==' ' || *ch=='\t') ch++;
1808 if (*ch=='\0' || *ch=='\n') goto err3;
1809 }
1810 if (valprs(ch,&v1)==NULL) goto err1;
1811 *vvs = v1;
1812 vvs++;
1813 }
1814 *vvs = -999.9;
1815 pfi->abvals[dim] = vals;
1816 pfi->grvals[dim] = vals;
1817 pfi->ab2gr[dim] = lev2gr;
1818 pfi->gr2ab[dim] = gr2lev;
1819 pfi->linear[dim] = 0;
1820 return (0);
1821
1822 err1:
1823 gaprnt (0,"Open Error: Invalid value in LEVELS data\n");
1824 free (vals);
1825 return (1);
1826
1827 err2:
1828 gaprnt (0,"Open Error: Unexpected EOF reading descriptor file\n");
1829 gaprnt (0," EOF occurred reading LEVELS values\n");
1830 free (vals);
1831 return (1);
1832
1833 err3:
1834 gaprnt (0,"Open Error: Blank Record found in LEVELS data\n");
1835 free (vals);
1836 return (1);
1837 }
1838
1839
1840 /* Process attribute metadata embedded in a comment record */
1841 /* Return codes: -1 is memory allocation error */
1842
ddfattr(char * ch,struct gafile * pfi)1843 int ddfattr (char *ch, struct gafile *pfi) {
1844
1845 int jj,len;
1846 char *varname,*attrname,*attrtype,attrvalue[512];
1847 struct gaattr *testattr,*attrib;
1848 varname = attrname = attrtype = NULL;
1849
1850 /* check for presence of attribute metadata */
1851 if ((ch=nxtwrd(ch))==NULL ) {
1852 gaprnt (2,"Warning: Attribute comment missing all required fields \n");
1853 return(0);
1854 }
1855
1856 /* get the variable name */
1857 len = 0;
1858 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
1859 varname = (char *)malloc(len+3);
1860 for (jj=0; jj<len; jj++) *(varname+jj) = *(ch+jj);
1861 *(varname+len) = '\0';
1862
1863 /* get the attribute type */
1864 if ( (ch=nxtwrd(ch))==NULL ) {
1865 gaprnt (2,"Warning: Attribute comment missing attribute type, name, and value \n");
1866 if (varname) { free (varname); varname = NULL; }
1867 return(0);
1868 }
1869 len = 0;
1870 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
1871 attrtype = (char *)malloc(len+3);
1872 for (jj=0; jj<len; jj++) *(attrtype+jj) = *(ch+jj);
1873 *(attrtype+len) = '\0';
1874
1875 /* check if attribute type matches list of valid types */
1876 if ((strncmp(attrtype,"String",6) != 0) &&
1877 (strncmp(attrtype,"Str",3) != 0) &&
1878 (strncmp(attrtype,"Byte",4) != 0) &&
1879 (strncmp(attrtype,"Int16",5) != 0) &&
1880 (strncmp(attrtype,"UInt16",6) != 0) &&
1881 (strncmp(attrtype,"Int32",5) != 0) &&
1882 (strncmp(attrtype,"UInt32",6) != 0) &&
1883 (strncmp(attrtype,"Float32",7) != 0) &&
1884 (strncmp(attrtype,"Float64",7) != 0) &&
1885 (strncmp(attrtype,"Url",3) != 0)) {
1886 sprintf (pout,"Warning: Attribute type \"%s\" is not one of the OPeNDAP standards:\n",attrtype);
1887 gaprnt (2,pout);
1888 sprintf (pout," String, Str, Byte, Int16, UInt16, Int32, UInt32, Float32, Float64, Url\n");
1889 gaprnt (2,pout);
1890 }
1891
1892 /* get the attribute name */
1893 if ( (ch=nxtwrd(ch))==NULL ) {
1894 gaprnt (2,"Warning: Attribute comment missing attribute name and value \n");
1895 if (varname) { free (varname); varname = NULL; }
1896 if (attrtype) { free (attrtype); attrtype = NULL; }
1897 return(0);
1898 }
1899 len = 0;
1900 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
1901 attrname = (char *)malloc(len+3);
1902 for (jj=0; jj<len; jj++) *(attrname+jj) = *(ch+jj);
1903 *(attrname+len) = '\0';
1904
1905 /* get the attribute value */
1906 if ( (ch=nxtwrd(ch))==NULL ) {
1907 gaprnt (2,"Warning: Attribute comment missing attribute value \n");
1908 if (varname) { free (varname); varname = NULL; }
1909 if (attrtype) { free (attrtype); attrtype = NULL; }
1910 if (attrname) { free (attrname); attrname = NULL; }
1911 return(0);
1912 }
1913 getstr (attrvalue, ch, 512);
1914
1915 /* Everything parsed OK, so add a link to the list of attributes */
1916 if (pfi->attr) { /* chain already exists */
1917 /* advance to end of chain */
1918 attrib=pfi->attr;
1919 while (attrib->next) attrib=attrib->next;
1920 /* initialize new link */
1921 if ((attrib->next = (struct gaattr *) calloc(1,sizeof(struct gaattr))) != NULL) {
1922 attrib = attrib->next; /* advance to new end of chain */
1923 attrib->varname = varname;
1924 attrib->name = attrname;
1925 attrib->type = attrtype;
1926 strcpy(attrib->value,attrvalue);
1927 attrib->next = NULL;
1928 } else {
1929 gaprnt (0,"Error: malloc failed for attribute metadata\n");
1930 if (varname) { free (varname); varname = NULL; }
1931 if (attrtype) { free (attrtype); attrtype = NULL; }
1932 if (attrname) { free (attrname); attrname = NULL; }
1933 return(-1);
1934 }
1935 }
1936 else {
1937 /* initialize first link */
1938 pfi->attr = (struct gaattr *) calloc(1,sizeof(struct gaattr));
1939 if (pfi->attr != NULL) {
1940 pfi->attr->varname = varname;
1941 pfi->attr->name = attrname;
1942 pfi->attr->type = attrtype;
1943 strcpy(pfi->attr->value,attrvalue);
1944 pfi->attr->next = NULL;
1945 attrib = pfi->attr;
1946 } else {
1947 gaprnt (0,"Error: malloc failed for attribute metadata\n");
1948 if (varname) { free (varname); varname = NULL; }
1949 if (attrtype) { free (attrtype); attrtype = NULL; }
1950 if (attrname) { free (attrname); attrname = NULL; }
1951 return(-1);
1952 }
1953 }
1954 return(0);
1955 }
1956
1957
1958 /* Allocate and initialize a gafile structure */
1959
getpfi(void)1960 struct gafile *getpfi (void) {
1961 struct gafile *pfi;
1962 int i;
1963
1964 pfi = (struct gafile *)malloc(sizeof(struct gafile));
1965 if (pfi==NULL) return (NULL);
1966
1967 pfi->type = 1; /* Assume grid unless told otherwise */
1968 pfi->tlpflg = 0; /* Assume file not circular */
1969 pfi->bswap = 0; /* Assume no byte swapping needed */
1970 pfi->seqflg = 0; /* Assume direct access */
1971 pfi->yrflg = 0; /* Assume south to north */
1972 pfi->zrflg = 0; /* Assume bottom to top */
1973 pfi->idxflg = 0; /* Assume binary */
1974 pfi->ncflg = 0; /* Assume not netcdf */
1975 pfi->bufrflg = 0; /* Assume not bufr */
1976 pfi->ncid = -999; /* No netcdf file open */
1977 pfi->sdid = -999; /* No hdfsds file open */
1978 pfi->fhdr = 0; /* Assume no file header */
1979 pfi->thdr=0; /*mf assume no time header mf*/
1980 pfi->xyhdr=0; /*mf assume no xyheader mf*/
1981 pfi->fseq = -999; /* No sequence number assigned */
1982 pfi->dhandle = -999; /* Assume not a gadods stn data set */
1983 pfi->packflg = 0; /* Assume data are not packed */
1984 pfi->undefattrflg = 0; /* Assume no undef attribute name given */
1985
1986 pfi->bufrinfo = NULL;
1987 pfi->bufrdset = NULL;
1988 pfi->attr = NULL;
1989 pfi->scattr = NULL;
1990 pfi->ofattr = NULL;
1991 pfi->undefattr = NULL;
1992 pfi->tempname = NULL;
1993 pfi->mnam = NULL;
1994 pfi->infile = NULL;
1995 pfi->rbuf = NULL;
1996 pfi->pbuf = NULL;
1997 pfi->bbuf = NULL;
1998 pfi->tstrt = NULL;
1999 pfi->tcnt = NULL;
2000 pfi->mfile = NULL;
2001 pfi->pvar1 = NULL;
2002 pfi->pindx = NULL;
2003 pfi->fnums = NULL;
2004 pfi->pchsub1 = NULL;
2005
2006 pfi->wrap = 0; /* Assume no wrapping */
2007 for (i=0; i<4; i++) pfi->dimoff[i] = 0;
2008 pfi->title[0] = '\0';
2009 pfi->grvals[0] = NULL;
2010 pfi->grvals[1] = NULL;
2011 pfi->grvals[2] = NULL;
2012 pfi->grvals[3] = NULL;
2013 pfi->grbgrd = -999;
2014 pfi->tmplat = 0;
2015 pfi->errcnt = 0;
2016 pfi->errflg = 0;
2017 pfi->ppflag = 0; /* Assume lat-lon grid */
2018 pfi->ppwrot = 0; /* Assume no wind rotataion */
2019 for (i=0; i<9; i++) pfi->ppi[i] = NULL;
2020 for (i=0; i<9; i++) pfi->ppf[i] = NULL;
2021 pfi->ppw = NULL;
2022 pfi->calendar=0;
2023 pfi->cray_ieee=0;
2024
2025 #if USESDF == 1
2026 pfi->sdf_ptr = NULL ;
2027 pfi->is_a_SDF = 0 ;
2028 pfi->first_sdf_ptr = NULL ;
2029 #endif
2030
2031 return (pfi);
2032 }
2033
2034 /* Free a gafile structure and associated storage. If the flag is
2035 true, DO NOT free the storage related to scaling transforms,
2036 since someone, somewhere, may still be pointing to that. */
2037
frepfi(struct gafile * pfi,int flag)2038 void frepfi (struct gafile *pfi, int flag) {
2039 struct gaattr *attrib, *nextattrib;
2040 struct gaindx *pindx;
2041 struct gachsub *pchsub,*pch2;
2042 int i;
2043 #if USESDF == 1
2044 void free_io_std(IO_STD**) ;
2045 #endif
2046
2047 if (pfi->pindx) {
2048 pindx = pfi->pindx;
2049 if (pindx->hipnt) free(pindx->hipnt);
2050 if (pindx->hfpnt) free(pindx->hfpnt);
2051 if (pindx->intpnt) free(pindx->intpnt);
2052 if (pindx->fltpnt) free(pindx->fltpnt);
2053 free (pindx);
2054 }
2055 pchsub = pfi->pchsub1;
2056 while (pchsub) {
2057 if (pchsub->ch) free(pchsub->ch);
2058 pch2 = pchsub->forw;
2059 free (pchsub);
2060 pchsub = pch2;
2061 }
2062 if (pfi->fnums) free(pfi->fnums);
2063 if (pfi->tstrt) free(pfi->tstrt);
2064 if (pfi->tcnt) free(pfi->tcnt);
2065 if (pfi->pvar1) free(pfi->pvar1);
2066 if (pfi->rbuf) free(pfi->rbuf);
2067 if (pfi->pbuf) free(pfi->pbuf);
2068 if (pfi->mnam) free(pfi->mnam);
2069 if (pfi->scattr) free(pfi->scattr);
2070 if (pfi->ofattr) free(pfi->ofattr);
2071 if (pfi->undefattr) free(pfi->undefattr);
2072 if (pfi->bufrinfo) free(pfi->bufrinfo);
2073 #ifndef STNDALN
2074 if (pfi->bufrdset) gabufr_close(pfi->bufrdset);
2075 #endif
2076 for (i=0; i<9; i++) if (pfi->ppi[i]) free(pfi->ppi[i]);
2077 for (i=0; i<9; i++) if (pfi->ppf[i]) free(pfi->ppf[i]);
2078 if (pfi->ppw) free(pfi->ppw);
2079 if (!flag) {
2080 for (i=0; i<4; i++) if (pfi->grvals[i]) free(pfi->grvals[i]);
2081 }
2082 if (pfi->attr) {
2083 for (attrib = pfi->attr; attrib != NULL; attrib = nextattrib) {
2084 nextattrib = attrib->next;
2085 if (attrib->varname) free(attrib->varname);
2086 if (attrib->type) free(attrib->type);
2087 if (attrib->name) free(attrib->name);
2088 free(attrib);
2089 }
2090 }
2091 #if USESDF == 1
2092 if (pfi->sdf_ptr) {
2093 if ((pfi->sdf_ptr) == (pfi->first_sdf_ptr)) {
2094 pfi->first_sdf_ptr = (IO_STD *) 0 ;
2095 }
2096 free_io_std(&(pfi->sdf_ptr)) ;
2097 }
2098 #endif
2099 free (pfi);
2100 }
2101
2102
2103 /* Routine to calculate or input the interpolation constants needed for
2104 the implicit interpolation from pre-projected grids to lat-lon. */
2105
gappcn(struct gafile * pfi,int pdefop1,int pdefop2)2106 int gappcn (struct gafile *pfi, int pdefop1, int pdefop2) {
2107 int size,i,j,ii,jj;
2108 float lat,lon,rii,rjj;
2109 float *dx, *dy, *dw, dum;
2110 float pi;
2111 int *ioff, rdw, rc, pnum, wflg;
2112
2113 dw=NULL;
2114 size = pfi->dnum[0]*pfi->dnum[1];
2115
2116 /* Allocate space needed for the ppi and ppf grids */
2117 if (pfi->ppflag != 8) {
2118 ;
2119 if ((pfi->ppi[0] = (int *)malloc(sizeof(int )*size)) == NULL) goto merr;
2120 if ((pfi->ppf[0] = (float *)malloc(sizeof(float)*size)) == NULL) goto merr;
2121 if ((pfi->ppf[1] = (float *)malloc(sizeof(float)*size)) == NULL) goto merr;
2122 if (pfi->ppwrot) {
2123 if ((pfi->ppw = (float *)malloc(sizeof(float)*size)) == NULL) goto merr;
2124 }
2125 }
2126
2127 /* pdef bilin */
2128 if (pfi->ppflag==7) {
2129 /* allocate memory for ppi values */
2130 if ((pfi->ppi[1] = (int *)malloc(sizeof(int)*size)) == NULL) goto merr;
2131
2132 if (pdefop1==2) { /* sequential -- read the 4-byte header */
2133 rc = fread(&rdw, sizeof(int), 1, pdfi);
2134 if (rc!=1) goto merr2;
2135 }
2136 /* read the grid of pdef ivals */
2137 rc = fread(pfi->ppf[0], sizeof(float), size, pdfi); if (rc!=size) goto merr2;
2138
2139 if (pdefop1==2) { /* sequential -- read the 4-byte footer and next header */
2140 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2141 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2142 }
2143 /* read the grid of pdef jvals */
2144 rc = fread(pfi->ppf[1], sizeof(float), size, pdfi); if (rc!=size) goto merr2;
2145
2146 if (pdefop1==2) { /* sequential -- read the 4-byte footer and next header */
2147 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2148 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2149 }
2150 /* read the grid of wind rotation vals */
2151 rc = fread(pfi->ppw, sizeof(float), size, pdfi); if (rc!=size) goto merr2;
2152
2153 /* byte swap if necessary */
2154 if ((pdefop2==2 && !BYTEORDER) ||
2155 (pdefop2==3 && BYTEORDER)) {
2156 gabswp (pfi->ppf[0],size);
2157 gabswp (pfi->ppf[1],size);
2158 gabswp (pfi->ppw,size);
2159 }
2160
2161 /* Fill grids of file offsets and weights (dx,dy) for pdef grid interpolation */
2162 ioff = pfi->ppi[0];
2163 dx = pfi->ppf[0];
2164 dy = pfi->ppf[1];
2165 dw = pfi->ppw;
2166 wflg = 0;
2167 for (j=0; j<pfi->dnum[1]; j++) {
2168 for (i=0; i<pfi->dnum[0]; i++) {
2169 if (*dx < 0.0) *ioff = -1;
2170 else {
2171 ii = (int)(*dx);
2172 jj = (int)(*dy);
2173 *dx = *dx - (float)ii;
2174 *dy = *dy - (float)jj;
2175 if (ii<1 || ii>pfi->ppisiz-1 || jj<1 || jj>pfi->ppjsiz-1) {
2176 *ioff = -1;
2177 } else {
2178 *ioff = (jj-1)*pfi->ppisiz + ii - 1;
2179 }
2180 }
2181 if (fabs(*dw) > 0.00001) wflg = 1;
2182 ioff++; dx++; dy++, dw++;
2183 }
2184 }
2185 pfi->ppwrot = wflg;
2186
2187 /* When pdef is a file, read in the offsets of the points
2188 to use and their weights, as well as the array of
2189 wind rotation values to use */
2190
2191 } else if (pfi->ppflag==8) {
2192 pnum = (int)(pfi->ppvals[0]+0.1);
2193 for (i=0; i<pnum; i++) {
2194 /* allocate memory and then read in the offsets */
2195 if ((pfi->ppi[i] = (int *)malloc(sizeof(int)*size)) == NULL) goto merr;
2196 if (pdefop1==2) { /* sequential -- header */
2197 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2198 }
2199 rc = fread(pfi->ppi[i], sizeof(int), size, pdfi); if (rc!=size) goto merr2;
2200 if (pdefop1==2) { /* sequential -- footer */
2201 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2202 }
2203
2204 /* allocate memory and then read in the weights */
2205 if ((pfi->ppf[i] = (float *)malloc(sizeof(int)*size)) == NULL) goto merr;
2206 if (pdefop1==2) { /* sequential -- header */
2207 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2208 }
2209 rc = fread(pfi->ppf[i], sizeof(float), size, pdfi); if (rc!=size) goto merr2;
2210 if (pdefop1==2) { /* sequential -- footer */
2211 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2212 }
2213
2214 /* byte swap if necessary */
2215 if ((pdefop2==2 && !BYTEORDER) ||
2216 (pdefop2==3 && BYTEORDER)) {
2217 gabswp ((float *)(pfi->ppi[i]),size);
2218 gabswp (pfi->ppf[i],size);
2219 }
2220 }
2221
2222 /* allocate memory and read in the wind rotation values */
2223 if ((pfi->ppw = (float *)malloc(sizeof(float)*size)) == NULL) goto merr;
2224 if (pdefop1==2) { /* sequential -- header */
2225 rc = fread(&rdw, sizeof(int), 1, pdfi); if (rc!=1) goto merr2;
2226 }
2227 rc = fread(pfi->ppw, sizeof(float), size, pdfi); if (rc!=size) goto merr2;
2228
2229 /* byte swap if necessary */
2230 if ((pdefop2==2 && !BYTEORDER) ||
2231 (pdefop2==3 && BYTEORDER)) {
2232 gabswp (pfi->ppw,size);
2233 }
2234
2235 /* set wind rotation flag */
2236 dw = pfi->ppw;
2237 wflg = 0;
2238 for (i=0; i<size; i++) {
2239 if (fabs(*dw) > 0.00001) wflg = 1;
2240 dw++;
2241 }
2242 pfi->ppwrot = wflg;
2243
2244 /* When a supported projection is specified,
2245 we calculate three constants at each lat-lon grid point:
2246 offset of the ij gridpoint, and the delta x and delta y
2247 values. */
2248
2249 } else {
2250 pi = acos(-1.0);
2251 ioff = pfi->ppi[0];
2252 dx = pfi->ppf[0];
2253 dy = pfi->ppf[1];
2254 if ( pfi->ppwrot) dw = pfi->ppw;
2255
2256 for (j=0; j<pfi->dnum[1]; j++) {
2257 lat = pfi->gr2ab[1](pfi->grvals[1],(float)(j+1));
2258 for (i=0; i<pfi->dnum[0]; i++) {
2259 lon = pfi->gr2ab[0](pfi->grvals[0],(float)(i+1));
2260 /* get i,j values in preprojected grid for each lat/lon point */
2261 if (pfi->ppflag==3) {
2262 if(pfi->ppwrot) { /* PDEF lccr */
2263 ll2lc (pfi->ppvals, lat, lon, &rii, &rjj, dw);
2264 } else { /* PDEF lcc */
2265 ll2lc (pfi->ppvals, lat, lon, &rii, &rjj, &dum);
2266 }
2267 } else if (pfi->ppflag==4) { /* PDEF eta.u */
2268 ll2eg (pfi->ppisiz,pfi->ppjsiz,pfi->ppvals, lon, lat, &rii, &rjj, dw);
2269 } else if (pfi->ppflag==5) { /* PDEF pse */
2270 ll2pse (pfi->ppisiz,pfi->ppjsiz,pfi->ppvals, lon, lat, &rii, &rjj);
2271 } else if (pfi->ppflag==6) { /* PDEF ops */
2272 ll2ops (pfi->ppvals, lon, lat, &rii, &rjj);
2273 } else { /* PDEF nps and sps */
2274 w3fb04(lat, -1.0*lon, pfi->ppvals[3], -1.0*pfi->ppvals[2], &rii, &rjj);
2275 rii = rii + pfi->ppvals[0]; /* Normalize based on pole point */
2276 rjj = rjj + pfi->ppvals[1];
2277 *dw = (pfi->ppvals[2]-lon) * pi/180.0; /* wind rotation amount */
2278 if (pfi->ppflag==2) *dw = pi - *dw;
2279 }
2280 ii = (int)rii;
2281 jj = (int)rjj;
2282 *dx = rii - (float)ii;
2283 *dy = rjj - (float)jj;
2284 if (ii<1 || ii>pfi->ppisiz-1 || jj<1 || jj>pfi->ppjsiz-1) {
2285 *ioff = -1;
2286 } else {
2287 *ioff = (jj-1)*pfi->ppisiz + ii - 1;
2288 }
2289 ioff++; dx++; dy++;
2290 if (pfi->ppwrot) dw++;
2291 }
2292 }
2293 }
2294 return(0);
2295
2296 merr:
2297 gaprnt (0,"Open Error: Memory allocation error in pdef handler\n");
2298 return(1);
2299
2300 merr2:
2301 gaprnt (0,"Open Error: I/O Error on pdef file read\n");
2302 return(1);
2303
2304 }
2305
w3fb04(float alat,float along,float xmeshl,float orient,float * xi,float * xj)2306 void w3fb04 (float alat, float along, float xmeshl, float orient,
2307 float *xi, float *xj) {
2308
2309 /*
2310 C
2311 C SUBPROGRAM: W3FB04 LATITUDE, LONGITUDE TO GRID COORDINATES
2312 C AUTHOR: MCDONELL,J. ORG: W345 DATE: 90-06-04
2313 C
2314 C ABSTRACT: CONVERTS THE COORDINATES OF A LOCATION ON EARTH FROM THE
2315 C NATURAL COORDINATE SYSTEM OF LATITUDE/LONGITUDE TO THE GRID (I,J)
2316 C COORDINATE SYSTEM OVERLAID ON A POLAR STEREOGRAPHIC MAP PRO-
2317 C JECTION TRUE AT 60 DEGREES N OR S LATITUDE. W3FB04 IS THE REVERSE
2318 C OF W3FB05.
2319 C
2320 C PROGRAM HISTORY LOG:
2321 C 77-05-01 J. MCDONELL
2322 C 89-01-10 R.E.JONES CONVERT TO MICROSOFT FORTRAN 4.1
2323 C 90-06-04 R.E.JONES CONVERT TO SUN FORTRAN 1.3
2324 C 93-01-26 B. Doty converted to C
2325 C
2326 C USAGE: CALL W3FB04 (ALAT, ALONG, XMESHL, ORIENT, XI, XJ)
2327 C
2328 C INPUT VARIABLES:
2329 C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES
2330 C ------ --------- -----------------------------------------------
2331 C ALAT ARG LIST LATITUDE IN DEGREES (<0 IF SH)
2332 C ALONG ARG LIST WEST LONGITUDE IN DEGREES
2333 C XMESHL ARG LIST MESH LENGTH OF GRID IN KM AT 60 DEG LAT(<0 IF SH)
2334 C (190.5 LFM GRID, 381.0 NH PE GRID,-381.0 SH PE GRID)
2335 C ORIENT ARG LIST ORIENTATION WEST LONGITUDE OF THE GRID
2336 C (105.0 LFM GRID, 80.0 NH PE GRID, 260.0 SH PE GRID)
2337 C
2338 C OUTPUT VARIABLES:
2339 C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES
2340 C ------ --------- -----------------------------------------------
2341 C XI ARG LIST I OF THE POINT RELATIVE TO NORTH OR SOUTH POLE
2342 C XJ ARG LIST J OF THE POINT RELATIVE TO NORTH OR SOUTH POLE
2343 C
2344 C SUBPROGRAMS CALLED:
2345 C NAMES LIBRARY
2346 C ------------------------------------------------------- --------
2347 C COS SIN SYSLIB
2348 C
2349 C REMARKS: ALL PARAMETERS IN THE CALLING STATEMENT MUST BE
2350 C REAL. THE RANGE OF ALLOWABLE LATITUDES IS FROM A POLE TO
2351 C 30 DEGREES INTO THE OPPOSITE HEMISPHERE.
2352 C THE GRID USED IN THIS SUBROUTINE HAS ITS ORIGIN (I=0,J=0)
2353 C AT THE POLE IN EITHER HEMISPHERE, SO IF THE USER'S GRID HAS ITS
2354 C ORIGIN AT A POINT OTHER THAN THE POLE, A TRANSLATION IS NEEDED
2355 C TO GET I AND J. THE GRIDLINES OF I=CONSTANT ARE PARALLEL TO A
2356 C LONGITUDE DESIGNATED BY THE USER. THE EARTH'S RADIUS IS TAKEN
2357 C TO BE 6371.2 KM.
2358 C
2359 C ATTRIBUTES:
2360 C LANGUAGE: SUN FORTRAN 1.4
2361 C MACHINE: SUN SPARCSTATION 1+
2362 C*/
2363
2364 static float d2r = 3.14159/180.0;
2365 static float earthr = 6371.2;
2366
2367 float re,xlat,wlong,r;
2368
2369 re = (earthr * 1.86603) / xmeshl;
2370 xlat = alat * d2r;
2371
2372 if (xmeshl>0.0) {
2373 wlong = (along + 180.0 - orient) * d2r;
2374 r = (re * cos(xlat)) / (1.0 + sin(xlat));
2375 *xi = r * sin(wlong);
2376 *xj = r * cos(wlong);
2377
2378 } else {
2379
2380 re = -re;
2381 xlat = -xlat;
2382 wlong = (along - orient) * d2r;
2383 r = (re * cos(xlat)) / (1.0+ sin(xlat));
2384 *xi = r * sin(wlong);
2385 *xj = -r * cos(wlong);
2386 }
2387 }
2388
2389 /* Lambert conformal conversion */
2390
ll2lc(float * vals,float grdlat,float grdlon,float * grdi,float * grdj,float * wrot)2391 void ll2lc (float *vals, float grdlat, float grdlon, float *grdi, float *grdj, float *wrot) {
2392
2393 /* Subroutine to convert from lat-lon to lambert conformal i,j.
2394 Provided by NRL Monterey; converted to C 6/15/94.
2395
2396 c SUBROUTINE: ll2lc
2397 c
2398 c PURPOSE: To compute i- and j-coordinates of a specified
2399 c grid given the latitude and longitude points.
2400 c All latitudes in this routine start
2401 c with -90.0 at the south pole and increase
2402 c northward to +90.0 at the north pole. The
2403 c longitudes start with 0.0 at the Greenwich
2404 c meridian and increase to the east, so that
2405 c 90.0 refers to 90.0E, 180.0 is the inter-
2406 c national dateline and 270.0 is 90.0W.
2407 c
2408 c INPUT VARIABLES:
2409 c
2410 c vals+0 latref: latitude at reference point (iref,jref)
2411 c
2412 c vals+1 lonref: longitude at reference point (iref,jref)
2413 c
2414 c vals+2 iref: i-coordinate value of reference point
2415 c
2416 c vals+3 jref: j-coordinate value of reference point
2417 c
2418 c vals+4 stdlt1: standard latitude of grid (S True lat)
2419 c
2420 c vals+5 stdlt2: second standard latitude of grid (only required
2421 c if igrid = 2, lambert conformal) (N True lat)
2422 c
2423 c vals+6 stdlon: standard longitude of grid (longitude that
2424 c points to the north)
2425 c
2426 c vals+7 delx: grid spacing of grid in x-direction
2427 c for igrid = 1,2,3 or 4, delx must be in meters
2428 c for igrid = 5, delx must be in degrees
2429 c
2430 c vals+8 dely: grid spacing (in meters) of grid in y-direction
2431 c for igrid = 1,2,3 or 4, delx must be in meters
2432 c for igrid = 5, dely must be in degrees
2433 c
2434 c grdlat: latitude of point (grdi,grdj)
2435 c
2436 c grdlon: longitude of point (grdi,grdj)
2437 c
2438 c grdi: i-coordinate(s) that this routine will generate
2439 c information for
2440 c
2441 c grdj: j-coordinate(s) that this routine will generate
2442 c information for
2443 c
2444 */
2445
2446 float pi, pi2, pi4, d2r, r2d, radius, omega4;
2447 float gcon,ogcon,H,deg,cn1,cn2,cn3,cn4,rih,xih,yih,rrih,check;
2448 float alnfix,alon,x,y,windrot;
2449 float latref,lonref,iref,jref,stdlt1,stdlt2,stdlon,delx,dely;
2450
2451 pi = 4.0*atan(1.0);
2452 pi2 = pi/2.0;
2453 pi4 = pi/4.0;
2454 d2r = pi/180.0;
2455 r2d = 180.0/pi;
2456 radius = 6371229.0;
2457 omega4 = 4.0*pi/86400.0;
2458
2459 latref = *(vals+0);
2460 lonref = *(vals+1);
2461 iref = *(vals+2);
2462 jref = *(vals+3);
2463 stdlt1 = *(vals+4);
2464 stdlt2 = *(vals+5);
2465 stdlon = *(vals+6);
2466 delx = *(vals+7);
2467 dely = *(vals+8);
2468
2469 /* case where standard lats are the same */
2470 /* corrected by Dan Geiszler of NRL; fabs of the
2471 lats was required for shem cases */
2472
2473 if(stdlt1 == stdlt2) {
2474 gcon = sin(d2r*(fabs(stdlt1)));
2475 } else {
2476 gcon = (log(sin((90.0-fabs(stdlt1))*d2r))
2477 -log(sin((90.0-fabs(stdlt2))*d2r)))
2478 /(log(tan((90.0-fabs(stdlt1))*0.5*d2r))
2479 -log(tan((90.0-fabs(stdlt2))*0.5*d2r)));
2480 }
2481 ogcon = 1.0/gcon;
2482 H = fabs(stdlt1)/(stdlt1); /* 1 for NHem, -1 for SHem */
2483 cn1 = sin((90.0-fabs(stdlt1))*d2r);
2484 cn2 = radius*cn1*ogcon;
2485 deg = (90.0-fabs(stdlt1))*d2r*0.5;
2486 cn3 = tan(deg);
2487 deg = (90.0-fabs(latref))*d2r*0.5;
2488 cn4 = tan(deg);
2489 rih = cn2*pow((cn4/cn3),gcon);
2490
2491 xih = rih*sin((lonref-stdlon)*d2r*gcon);
2492 yih = -rih*cos((lonref-stdlon)*d2r*gcon)*H;
2493 deg = (90.0-grdlat*H)*0.5*d2r;
2494 cn4 = tan(deg);
2495 rrih = cn2*pow((cn4/cn3),gcon);
2496 check = 180.0-stdlon;
2497 alnfix = stdlon+check;
2498 alon = grdlon+check;
2499
2500 while (alon< 0.0) alon = alon+360.0;
2501 while (alon>360.0) alon = alon-360.0;
2502
2503 deg = (alon-alnfix)*gcon*d2r;
2504 x = rrih*sin(deg);
2505 y = -rrih*cos(deg)*H;
2506 *grdi = iref + (x-xih)/delx;
2507 *grdj = jref + (y-yih)/dely;
2508 /* mf 20040630 -- use ftp://grads.iges.org/grads/src/grib212.f to calc rotation factor */
2509 windrot=gcon*(stdlon-grdlon)*d2r;
2510 *wrot=windrot;
2511 }
2512
2513 /* NMC eta ll to xy map */
2514
ll2eg(int im,int jm,float * vals,float grdlon,float grdlat,float * grdi,float * grdj,float * alpha)2515 void ll2eg (int im, int jm, float *vals, float grdlon, float grdlat,
2516 float *grdi, float *grdj, float *alpha) {
2517
2518 /* Subroutine to convert from lat-lon to NMC eta i,j.
2519
2520 Provided by Eric Rogers NMC;
2521 Converted to C 3/29/95 by Mike Fiorino
2522 Modified 9/2004 by J.M.Adams to correct grdi/j and alpha calculations
2523
2524 c SUBROUTINE: ll2eg
2525 c
2526 c PURPOSE: To compute i- and j-coordinates of a specified
2527 c grid given the latitude and longitude points.
2528 c All latitudes in this routine start
2529 c with -90.0 at the south pole and increase
2530 c northward to +90.0 at the north pole. The
2531 c longitudes start with 0.0 at the Greenwich
2532 c meridian and increase to the east, so that
2533 c 90.0 refers to 90.0E, 180.0 is the inter-
2534 c national dateline and 270.0 is 90.0W.
2535 c
2536 c INPUT VARIABLES:
2537 c
2538 c vals+0 tlm0d: longitude of the reference center point
2539 c vals+1 tph0d: latitude of the reference center point
2540 c vals+2 dlam: dlon grid increment in deg
2541 c vals+3 dphi: dlat grid increment in deg
2542 c
2543 c grdlat: latitude of point (grdi,grdj)
2544 c grdlon: longitude of point (grdi,grdj)
2545 c grdi: i-coordinate(s) that this routine will generate
2546 c information for
2547 c grdj: j-coordinate(s) that this routine will generate
2548 c information for
2549 */
2550
2551 float pi,d2r,r2d, earthr;
2552 float tlm0d,tph0d,dlam,dphi;
2553 float phi,lam,lam0,phi0;
2554 float x,y,z,xx,bigphi,biglam;
2555 float dlmd,dphd,wbd,sbd;
2556
2557 pi = 3.141592654;
2558 d2r = pi/180.0;
2559 r2d = 1.0/d2r;
2560 earthr = 6371.2;
2561
2562 tlm0d = -*(vals+0); /* convert + W to + E, the grads standard for longitude */
2563 tph0d = *(vals+1);
2564 dlam = (*(vals+2))*0.5;
2565 dphi = (*(vals+3))*0.5;
2566
2567 /* convert to radians */
2568 phi = grdlat*d2r; /* grid latitude */
2569 lam = -grdlon*d2r; /* grid longitude, convert +W to +E, the grads standard */
2570 phi0 = tph0d*d2r; /* center latitude */
2571 lam0 = tlm0d*d2r; /* center longitude */
2572
2573 /* Transform grid lat/lon */
2574 x = cos(phi0)*cos(phi)*cos(lam-lam0)+sin(phi0)*sin(phi);
2575 y = -cos(phi)*sin(lam-lam0);
2576 z = -sin(phi0)*cos(phi)*cos(lam-lam0)+cos(phi0)*sin(phi);
2577 biglam = atan2(y,x)/d2r; /* transformed lon in degrees */
2578 bigphi = atan2(z,sqrt(x*x+y*y))/d2r; /* transformed lat in degrees */
2579
2580 /* Convert transformed lat/lon -> i,j */
2581 dlmd = (*(vals+2));
2582 dphd = (*(vals+3));
2583 wbd = (-1)*0.5*(im-1)*dlmd; /* western boundary of transformed grid */
2584 sbd = (-1)*0.5*(jm-1)*dphd; /* southern boundary of transformed grid */
2585 *grdi = 1.0 + (biglam-wbd)/dlmd;
2586 *grdj = 1.0 + (bigphi-sbd)/dphd;
2587
2588 /* params for wind rotation alpha, alpha>0 ==> counter clockwise rotation */
2589 xx=sin(phi0)*sin(biglam*d2r)/cos(phi);
2590 if (xx < -1.0) xx = -1.0;
2591 else if (xx > 1.0) xx = 1.0;
2592 *alpha = (-1)*asin(xx);
2593
2594 }
2595
ll2pse(int im,int jm,float * vals,float lon,float lat,float * grdi,float * grdj)2596 void ll2pse (int im, int jm, float *vals, float lon, float lat,
2597 float *grdi, float *grdj) {
2598
2599
2600 /* Convert from geodetic latitude and longitude to polar stereographic
2601 grid coordinates. Follows mapll by V. J. Troisi. */
2602 /* Conventions include that slat and lat must be absolute values */
2603 /* The hemispheres are controlled by the sgn parameter */
2604 /* Bob Grumbine 15 April 1994. */
2605
2606 const float rearth = 6378.273e3;
2607 const float eccen2 = 0.006693883;
2608 const float pi = 3.141592654;
2609
2610 float cdr, alat, along, e, e2;
2611 float t, x, y, rho, sl, tc, mc;
2612 float slat,slon,xorig,yorig,sgn,polei,polej,dx,dy;
2613
2614 slat=*(vals+0);
2615 slon=*(vals+1);
2616 polei=*(vals+2);
2617 polej=*(vals+3);
2618 dx=*(vals+4)*1000;
2619 dy=*(vals+5)*1000;
2620 sgn=*(vals+6);
2621
2622 xorig = -polei*dx;
2623 yorig = -polej*dy;
2624
2625 cdr = 180./pi;
2626 alat = lat/cdr;
2627 along = lon/cdr;
2628 e2 = eccen2;
2629 e = sqrt(eccen2);
2630
2631 if ( fabs(lat) > 90.) {
2632 *grdi = -1;
2633 *grdj = -1;
2634 return;
2635 }
2636 else {
2637 t = tan(pi/4. - alat/2.) /
2638 pow( (1.-e*sin(alat))/(1.+e*sin(alat)) , e/2.);
2639
2640 if ( fabs(90. - slat) < 1.E-3) {
2641 rho = 2.*rearth*t/
2642 pow( pow(1.+e,1.+e) * pow(1.-e,1.-e) , e/2.);
2643 }
2644 else {
2645 sl = slat/cdr;
2646 tc = tan(pi/4.-sl/2.) /
2647 pow( (1.-e*sin(sl))/(1.+e*sin(sl)), (e/2.) );
2648 mc = cos(sl)/ sqrt(1.-e2*sin(sl)*sin(sl) );
2649 rho = rearth * mc*t/tc;
2650 }
2651
2652 x = rho*sgn*cos(sgn*(along+slon/cdr));
2653 y = rho*sgn*sin(sgn*(along+slon/cdr));
2654
2655 *grdi = (x - xorig)/dx+1;
2656 *grdj = (y - yorig)/dy+1;
2657
2658 return;
2659 }
2660
2661 }
2662
ll2ops(float * vals,float lni,float lti,float * grdi,float * grdj)2663 void ll2ops(float *vals, float lni, float lti, float *grdi, float *grdj) {
2664
2665 const float radius = 6371229.0 ;
2666
2667 float stdlat, stdlon, xref, yref, xiref, yjref, delx , dely;
2668 float plt,pln;
2669 double pi180,c1,c2,c3,c4,c5,c6,arg2a,bb,plt1,alpha, pln1,plt90,argu1,argu2;
2670 double hsign,glor,rstdlon,glolim,facpla,x,y;
2671
2672 stdlat = *(vals+0);
2673 stdlon = *(vals+1);
2674 xref = *(vals+2);
2675 yref = *(vals+3);
2676 xiref = *(vals+4);
2677 yjref = *(vals+5);
2678 delx = *(vals+6);
2679 dely = *(vals+7);
2680
2681 c1=1.0 ;
2682 pi180 = asin(c1)/90.0;
2683
2684 /* set flag for n/s hemisphere and convert longitude to <0 ; 360> interval */
2685 if(stdlat >= 0.0) {
2686 hsign= 1.0 ;
2687 } else {
2688 hsign=-1.0 ;
2689 }
2690
2691 /* set flag for n/s hemisphere and convert longitude to <0 ; 360> interval */
2692 glor=lni ;
2693 if(glor <= 0.0) glor=360.0+glor ;
2694 rstdlon=stdlon;
2695 if(rstdlon < 0.0) rstdlon=360.0+stdlon;
2696
2697 /* test for a n/s pole case */
2698 if(stdlat == 90.0) {
2699 plt=lti ;
2700 pln=fmod(glor+270.0,360.0) ;
2701 goto l2000;
2702 }
2703
2704 if(stdlat == -90.0) {
2705 plt=-lti ;
2706 pln=fmod(glor+270.0,360.0) ;
2707 goto l2000;
2708 }
2709
2710 /* test for longitude on 'greenwich or date line' */
2711 if(glor == rstdlon) {
2712 if(lti > stdlat) {
2713 plt=90.0-lti+stdlat;
2714 pln=90.0;
2715 } else {
2716 plt=90.0-stdlat+lti;
2717 pln=270.0;;
2718 }
2719 goto l2000;
2720 }
2721
2722 if(fmod(glor+180.0,360.0) == rstdlon) {
2723 plt=stdlat-90.0+lti;
2724 if(plt < -90.0) {
2725 plt=-180.0-plt;
2726 pln=270.0;
2727 } else {
2728 pln= 90.0;
2729 }
2730 goto l2000 ;
2731 }
2732
2733 /* determine longitude distance relative to rstdlon so it belongs to
2734 the absolute interval 0 - 180 */
2735 argu1 = glor-rstdlon;
2736 if(argu1 > 180.0) argu1 = argu1-360.0;
2737 if(argu1 < -180.0) argu1 = argu1+360.0;
2738
2739 /* 1. get the help circle bb and angle alpha (legalize arguments) */
2740
2741 c2=lti*pi180 ;
2742 c3=argu1*pi180 ;
2743 arg2a = cos(c2)*cos(c3) ;
2744 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = max1(arg2a,-c1) */
2745 if( c1 < arg2a ) arg2a = c1 ; /* min1(arg2a, c1) */
2746 bb = acos(arg2a) ;
2747
2748 c4=hsign*lti*pi180 ;
2749 arg2a = sin(c4)/sin(bb) ;
2750 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2751 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2752 alpha = asin(arg2a) ;
2753
2754 /* 2. get plt and pln (still legalizing arguments) */
2755 c5=stdlat*pi180 ;
2756 c6=hsign*stdlat*pi180 ;
2757 arg2a = cos(c5)*cos(bb) + sin(c6)*sin(c4) ;
2758 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2759 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2760 plt1 = asin(arg2a) ;
2761
2762 arg2a = sin(bb)*cos(alpha)/cos(plt1) ;
2763
2764 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2765 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2766 pln1 = asin(arg2a) ;
2767
2768
2769 /* test for passage of the 90 degree longitude (duallity in pln)
2770 get plt for which pln=90 when lti is the latitude */
2771 arg2a = sin(c4)/sin(c6) ;
2772 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2773 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2774 plt90 = asin(arg2a) ;
2775
2776 /* get help arc bb and angle alpha */
2777 arg2a = cos(c5)*sin(plt90) ;
2778 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2779 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2780 bb = acos(arg2a) ;
2781
2782 arg2a = sin(c4)/sin(bb) ;
2783 if( -c1 > arg2a ) arg2a = -c1 ; /* arg2a = dmax1(arg2a,-c1) */
2784 if( c1 < arg2a ) arg2a = c1 ; /* arg2a = dmin1(arg2a, c1) */
2785 alpha = asin(arg2a) ;
2786
2787 /* get glolim - it is nesc. to test for the existence of solution */
2788 argu2 = cos(c2)*cos(bb) / (1.-sin(c4)*sin(bb)*sin(alpha)) ;
2789 if( fabs(argu2) > c1 ) {
2790 glolim = 999.0;
2791 } else {
2792 glolim = acos(argu2)/pi180;
2793 }
2794
2795 /* modify (if nesc.) the pln solution */
2796 if( ( fabs(argu1) > glolim && lti <= stdlat ) || ( lti > stdlat ) ) {
2797 pln1 = pi180*180.0 - pln1;
2798 }
2799
2800 /* the solution is symmetric so the direction must be if'ed */
2801 if(argu1 < 0.0) {
2802 pln1 = -pln1;
2803 }
2804
2805 /* convert the radians to degrees */
2806 plt = plt1/pi180 ;
2807 pln = pln1/pi180 ;
2808
2809 /* to obtain a rotated value (ie so x-axis in pol.ste. points east)
2810 add 270 to longitude */
2811 pln=fmod(pln+270.0,360.0) ;
2812
2813 l2000:
2814
2815 /*
2816 c this program convert polar stereographic coordinates to x,y ditto
2817 c longitude: 0 - 360 ; positive to the east
2818 c latitude : -90 - 90 ; positive for northern hemisphere
2819 c it is assumed that the x-axis point towards the east and
2820 c corresponds to longitude = 0
2821 c
2822 c tsp 20/06-89
2823 c
2824 c constants and functions
2825 */
2826 facpla = radius*2.0/(1.0+sin(plt*pi180))*cos(plt*pi180);
2827 x = facpla*cos(pln*pi180) ;
2828 y = facpla*sin(pln*pi180) ;
2829
2830 *grdi=(x-xref)/delx + xiref;
2831 *grdj=(y-yref)/dely + yjref;
2832
2833 return;
2834
2835 }
2836
2837
2838 #if USESDF == 1
2839 #ifdef STNDALN
2840
2841 #define Success 1
2842 #define Failure 0
2843
2844 /* initialize a netCDF standard file structure */
2845 int
init_io_std(std_ptr)2846 init_io_std (std_ptr)
2847 IO_STD **std_ptr; /* pointer to netCDF data structure */
2848 {
2849 int init_dim_info ();
2850 void free_io_std ();
2851
2852 if ((*std_ptr) != NULL)
2853 free_io_std (std_ptr);
2854
2855 if (((*std_ptr) = (IO_STD *) malloc (sizeof (IO_STD))) == NULL)
2856 {
2857 printf ("Could not allocate memory for netCDF/HDF-SDS structure.\n");
2858 return Failure;
2859 }
2860
2861 /* initialize contents of structure */
2862 (*std_ptr)->cdfid = -1;
2863 (*std_ptr)->ndims = 0;
2864 (*std_ptr)->nvars = 0;
2865 (*std_ptr)->ngatts = 0;
2866 (*std_ptr)->recdim = -1;
2867 (*std_ptr)->time_type = CDC;
2868
2869 (*std_ptr)->first_gattr = NULL;
2870
2871 /* don't create variable list quite yet */
2872 (*std_ptr)->var = NULL;
2873
2874 /* intialize dimension information */
2875 if (init_dim_info ((*std_ptr)->dimids, (*std_ptr)->dimnam,
2876 (*std_ptr)->dimsiz) != Success)
2877 return Failure;
2878
2879 return Success;
2880 }
2881
2882 /* initialize dimension information */
2883 int
init_dim_info(dimids,dimnam,dimsiz)2884 init_dim_info (dimids, dimnam, dimsiz)
2885 int *dimids; /* array of dimension IDs */
2886 char dimnam[MAX_NC_DIMS][MAX_NC_NAME + 1]; /* array of dimension
2887 * names */
2888 long *dimsiz; /* array of dimension sizes */
2889 {
2890 int i;
2891
2892 /* set all dimension information to bogus values */
2893 for (i = 0; i < MAX_NC_DIMS; i++)
2894 {
2895 dimids[i] = -1;
2896 dimnam[i][0] = '\0';
2897 dimsiz[i] = -1;
2898 }
2899
2900 return Success;
2901 }
2902
2903 /* free a netCDF standard file structure */
2904 void
free_io_std(std_ptr)2905 free_io_std (std_ptr)
2906 IO_STD **std_ptr;
2907 {
2908 void free_var_info ();
2909 int free_netcdf_att_list ();
2910
2911 if (*std_ptr != NULL)
2912 {
2913 /* if the variable list exists, free it */
2914 if ((*std_ptr)->var != NULL)
2915 free_var_info (&((*std_ptr)->var));
2916
2917 /* free the global attributes if they have been allocated */
2918 if ((*std_ptr)->first_gattr != NULL)
2919 free_netcdf_att_list (&((*std_ptr)->first_gattr));
2920
2921 /* now free the IO_STD structure itself */
2922 free (*std_ptr);
2923 *std_ptr = NULL;
2924 }
2925
2926 return;
2927 }
2928
2929 /* free a variable structure */
2930 void
free_var_info(var)2931 free_var_info (var)
2932 VAR_INFO **var;
2933 {
2934 void free_var_info ();
2935 int free_netcdf_att_list ();
2936
2937 /* go through list to the end and free all variable structures */
2938 if ((*var)->next != NULL)
2939 free_var_info (&((*var)->next));
2940
2941 if ((*var)->first_vattr != NULL)
2942 free_netcdf_att_list (&((*var)->first_vattr));
2943
2944 /* if data space has been allocated, free it */
2945 if ((*var)->data != NULL)
2946 free ((*var)->data);
2947
2948 /* now free the variable structure itself */
2949 free (*var);
2950 *var = NULL;
2951
2952 return;
2953 }
2954
2955 /* Free a netCDF attribute list. */
2956 int
free_netcdf_att_list(first_attr)2957 free_netcdf_att_list (first_attr)
2958 struct attrib_list **first_attr; /* First attribute in list. */
2959 {
2960 struct attrib_list *temp_attr,
2961 *save_attr;
2962
2963 temp_attr = *first_attr;
2964 while (temp_attr != NULL)
2965 {
2966 save_attr = temp_attr->next;
2967 if (temp_attr->data != NULL)
2968 free (temp_attr->data);
2969 free (temp_attr);
2970 temp_attr = save_attr;
2971 }
2972 *first_attr = NULL;
2973 return Success;
2974 }
2975
2976 #endif
2977 #endif
2978