source: lmdz_wrf/trunk/WRFV3/external/io_grib1/grib1_routines.c @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 32.5 KB
Line 
1/*
2****************************************************************************
3*
4*  Routines for indexing, reading and writing grib files.  Routines
5*    are designed to be called by Fortran.
6*
7*  All routines return 0 for success, 1 for failure, unless otherwise noted.
8*
9*  Todd Hutchinson
10*  WSI
11*  05/17/2005
12*
13****************************************************************************
14*/
15
16#include "grib1_routines.h"
17#include "gridnav.h"
18#include "wrf_projection.h"
19#include <stdio.h>
20#include <stdlib.h>
21#include <ctype.h>
22#include <limits.h>
23
24char *trim (char *str);
25int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid);
26int index_times(GribInfo *gribinfo, Times *times);
27int find_time(Times *times, char valid_time[15]);
28int get_gridnav_projection(int wrf_projection);
29int get_byte(int input_int, int bytenum);
30
31/*
32 * Allocate space for the fileindex structure
33 */
34
35int ALLOC_INDEX_FILE(FileIndex *fileindex)
36{
37  int status = 0;
38 
39  fileindex->gribinfo = (GribInfo *)malloc(sizeof(GribInfo));
40  if (fileindex->gribinfo == NULL) {
41    fprintf(stderr,"Allocating fileindex->gribinfo failed.\n");
42    status = 1;
43    return status;
44  }
45 
46  fileindex->metadata = (MetaData *)malloc(sizeof(MetaData));
47  if (fileindex->metadata == NULL) {
48    fprintf(stderr,"Allocating fileindex->metadata failed.\n");
49    status = 1;
50    return status;
51  }
52  fileindex->metadata->elements = NULL;
53 
54  fileindex->times = (Times *)malloc(sizeof(Times));
55  if (fileindex->times == NULL) {
56    fprintf(stderr,"Allocating fileindex->times failed.\n");
57    status = 1;
58    return status;
59  }
60  fileindex->times->elements = NULL;
61 
62  return status;
63}
64
65
66void FREE_INDEX_FILE(FileIndex *fileindex)
67{
68  int status = 0;
69 
70  rg_free_gribinfo_elements(fileindex->gribinfo);
71  free(fileindex->gribinfo);
72
73  free(fileindex->metadata->elements);
74  free(fileindex->metadata);
75
76  free(fileindex->times->elements);
77  free(fileindex->times);
78
79}
80
81
82int INDEX_FILE(int *fid, FileIndex *fileindex)
83{
84
85  int status;
86  /* Index the grib records */
87
88  status = rg_setup_gribinfo_i(fileindex->gribinfo,*fid,1);
89  if (status < 0) {
90    fprintf(stderr,"Error setting up gribinfo structure.\n");
91    return 1;
92  }
93
94  /* Index the metadata section */
95
96  status = index_metadata(fileindex->gribinfo, fileindex->metadata, *fid);
97  if (status != 0) {
98    fprintf(stderr,"Error setting up metadata structure.\n");
99    return 1;
100  }
101
102  /* Setup a list of times based on times in grib records */
103
104  status = index_times(fileindex->gribinfo, fileindex->times);
105  if (status != 0) {
106    fprintf(stderr,"Error indexing times in grib file.\n");
107    return 1;
108  }
109 
110  return 0;
111}
112
113
114int GET_FILEINDEX_SIZE(int *size)
115{
116  *size = sizeof(FileIndex);
117  return *size;
118}
119
120
121int GET_NUM_TIMES(FileIndex *fileindex, int *numtimes)
122{
123  *numtimes = (fileindex->times)->num_elements;
124  return *numtimes;
125}
126
127
128int GET_TIME(FileIndex *fileindex, int *idx, char time[])
129{
130  int num_times;
131  int year, month, day, minute, hour, second;
132  char time2[100];
133
134  num_times = GET_NUM_TIMES(fileindex,&num_times);
135  if (*idx > num_times) 
136    {
137      fprintf(stderr,"Tried to get time %d, but only %d times exist\n",
138              *idx, num_times);
139      return 1;
140    }
141
142  strcpy(time,fileindex->times->elements[*idx-1].valid_time);
143 
144  /* Reformat time to meet WRF time format */
145
146  sscanf(time, "%4d%2d%2d%2d%2d%2d", 
147         &year, &month, &day, &hour, &minute, &second);
148  sprintf(time2, "%04d-%02d-%02d_%02d:%02d:%02d", 
149          year, month, day, hour, minute, second);
150  strncpy(time,time2,19);
151
152  return 0;
153}
154
155
156int GET_LEVEL1(FileIndex *fileindex, int *idx, int *level1)
157{
158
159  *level1 = (fileindex->gribinfo)->elements[*idx].usHeight1;
160
161  return *level1;
162
163}
164
165
166int GET_LEVEL2(FileIndex *fileindex, int *idx, int *level2)
167{
168
169  *level2 = (fileindex->gribinfo)->elements[*idx].usHeight2;
170
171  return *level2;
172
173}
174
175
176int index_metadata(GribInfo *gribinfo, MetaData *metadata, int fid)
177{
178  int status=0;
179  int end;
180  char string[11];
181  int found_metadata=0;
182  int idx;
183  int pos;
184  int fileend;
185  int seekpos;
186  int bytesread;
187  char line[1000];
188  char element[100],datestr[100],varname[100];
189  char value[1000];
190  int incomment;
191  int charidx;
192  int elemidx=0;
193  FILE *stream;
194
195
196  /* Associate a FILE *stream with the file id */
197  stream = fdopen(fid,"r");
198  if (stream == NULL) 
199    {
200      perror("Error associating stream with file descriptor");
201      status = -1;
202      return status;
203    }
204 
205  /*
206   * First, set the position to end of grib data (the
207   *   metadata section comes after the grib data).
208   */
209  idx = rg_num_elements(gribinfo) - 1;
210  end = rg_get_end(gribinfo,idx);
211  pos = end + 1;
212  fileend = lseek(fid,0,SEEK_END);
213
214  /*
215   * Now, start searching for metadata
216   */
217  while (pos < (fileend - 10))
218    {
219      seekpos = pos;
220      pos = lseek(fid,seekpos,SEEK_SET);
221      if (pos != seekpos) 
222        {
223          fprintf(stderr,"Error seeking %d bytes in file\n",end);
224          perror("");
225          return 1;
226        }
227
228      bytesread = read(fid,string,10);
229      if (bytesread != 10) 
230        {
231          fprintf(stderr,"Invalid read, pos: %d :\n",pos);
232          perror("");
233          pos += 1;
234          continue;
235        }
236
237      if (strncmp(string,"<METADATA>",10) == 0) 
238        {
239          /* We found it, so break out ! */
240          found_metadata = 1;
241          break;
242        }
243      pos += 1;
244    }
245
246
247  /* Now, read metadata, line by line */
248  incomment = 0;
249  while(fgets(line,1000,stream) != NULL)
250    {
251      trim(line);
252
253      /* Set comment flag, if we found a comment */
254      if (strncmp(line,"<!--",4) == 0)
255        {
256          incomment = 1;
257          strcpy(line,line+4);
258        }
259     
260      /* Search for end of comment */
261      if (incomment)
262        {
263          charidx = 0;
264          while (charidx < strlen(line)) 
265            {
266             
267              if (strncmp(line+charidx,"-->",3) == 0)
268                {
269                  strcpy(line,line+charidx+3);
270                  incomment = 0;
271                  break;
272                }
273              else
274                {
275                  charidx++;
276                }
277            }
278        }
279
280      if (incomment) continue;
281 
282
283      /* Check for end of metadata */
284      if (strncmp(line,"</METADATA>",11) == 0) 
285        {
286          /* We found end of data, so, break out */
287          break;
288        }
289
290      /* Skip blank lines */
291      if (strlen(line) == 0) continue;
292
293
294      /* Parse line */
295      trim(line);
296      strcpy(element,"none");
297      strcpy(datestr,"none");
298      strcpy(varname,"none");
299      strcpy(value,"none");
300      if (sscanf(line,"%[^;=];%[^;=];%[^;=]=%[^\n]",varname,datestr,
301                 element,value) == 4)
302        {
303        }
304      else if (sscanf(line,"%[^;=];%[^;=]=%[^\n]",datestr,element,value) == 3)
305        {
306          strcpy(varname,"none");
307        }
308      else if (sscanf(line,"%[^;=]=%[^\n]",element,value) == 2)
309        {
310          strcpy(varname,"none");
311          strcpy(datestr,"none");
312        }
313      else 
314        {
315          strcpy(varname,"none");
316          strcpy(datestr,"none");
317          strcpy(element,"none");
318          strcpy(value,"none");
319          fprintf(stderr,"Invalid line in metadata: \n%s",line);
320        }
321
322      trim(varname);
323      trim(datestr);
324      trim(element);
325      trim(value);
326
327      metadata->elements = 
328        (MetaData_Elements *)realloc( metadata->elements, 
329                                     (elemidx+1)*sizeof(MetaData_Elements) );
330      strcpy(metadata->elements[elemidx].VarName,varname);
331      strcpy(metadata->elements[elemidx].DateStr,datestr);
332      strcpy(metadata->elements[elemidx].Element,element);
333      strcpy(metadata->elements[elemidx].Value,value);
334
335      elemidx++;
336       
337    }
338
339  metadata->num_elements = elemidx;
340 
341  return 0;
342}
343
344
345
346
347int index_times(GribInfo *gribinfo, Times *times)
348{
349  int idx;
350  int status;
351  int numtimes=0;
352  int date;
353  char valid_time[15];
354  char tmp[15];
355  int swapped;
356
357  times->num_elements = 0;
358
359  /* Loop through elements, and build list of times */
360
361  for (idx=0; idx < gribinfo->num_elements; idx++) 
362    {
363      /* Calculate valid time */
364      status = rg_get_valid_time(gribinfo,idx,valid_time);
365      if (status != 0) 
366        {
367          fprintf(stderr,"Could not retrieve valid time for index: %d\n",idx);
368          continue;
369        }
370
371      /*
372       * Check if this time is already contained in times
373       *  If not, allocate space for it, and add it to list
374       */
375      if (find_time(times,valid_time) < 0) 
376        {
377          times->num_elements++;
378          times->elements = 
379            (Times_Elements *)
380            realloc(times->elements,times->num_elements*sizeof(Times_Elements));
381          if (times->elements == NULL) 
382            {
383              fprintf(stderr,"Allocating times->elements failed.\n");
384              status = 1;
385              return status;
386            }
387          strcpy(times->elements[times->num_elements - 1].valid_time,valid_time);
388        }
389    }
390
391  /* Sort times */
392  swapped = 1;
393  while (swapped)
394    {
395      swapped=0;
396      for (idx=1; idx < times->num_elements; idx++) 
397        {
398          if (strcmp(times->elements[idx-1].valid_time,
399                     times->elements[idx].valid_time) > 0)
400            {
401              strcpy(tmp,times->elements[idx-1].valid_time);
402              strcpy(times->elements[idx-1].valid_time,
403                     times->elements[idx].valid_time);
404              strcpy(times->elements[idx].valid_time, tmp);
405              swapped = 1;
406            }
407        }
408    }
409
410  return 0;
411}
412
413
414
415int find_time(Times *times, char valid_time[15])
416{
417  int idx;
418  int found_elem = -1;
419
420  for (idx = 0; idx < times->num_elements; idx++) 
421    {
422      if (strcmp(times->elements[idx].valid_time,valid_time) == 0) 
423        {
424          found_elem = idx;
425          break;
426        }
427    }
428
429  return found_elem;
430
431}
432
433
434int GET_METADATA_VALUE(FileIndex *fileindex, char ElementIn[], 
435                       char DateStrIn[], char VarNameIn[], char Value[], 
436                       int *stat, int strlen1, int strlen2, int strlen3, 
437                       int strlen4, int strlen5)
438{
439  int elemidx;
440  int elemnum;
441  char VarName[200];
442  char DateStr[200];
443  char Element[200];
444  int Value_Len;
445
446  *stat = 0;
447
448  strncpy(Element,ElementIn,strlen2);
449  Element[strlen2] = '\0';
450  strncpy(DateStr,DateStrIn,strlen3);
451  DateStr[strlen3] = '\0';
452  strncpy(VarName,VarNameIn,strlen4);
453  VarName[strlen4] = '\0';
454  Value_Len = strlen5;
455
456  elemnum = -1;
457  for (elemidx = 0; elemidx < fileindex->metadata->num_elements; elemidx++)
458    {
459      if (strcmp(Element,fileindex->metadata->elements[elemidx].Element) == 0)
460        {
461          if (strcmp(DateStr,fileindex->metadata->elements[elemidx].DateStr) 
462              == 0)
463            {
464              if (strcmp(VarName,
465                         fileindex->metadata->elements[elemidx].VarName) == 0)
466                {
467                  elemnum = elemidx;
468                  break;
469                }
470            }
471        }
472    }
473
474  if (elemnum != -1)
475    {
476      strncpy(Value, fileindex->metadata->elements[elemnum].Value, Value_Len);
477    } 
478  else
479    {
480      strncpy(Value, "none", Value_Len);
481      *stat = 1;
482    }
483
484  /*
485   * Pad end of string with one space.  This allows Fortran internal
486   *   read function to work properly.
487   */
488  if (strlen(Value) < Value_Len) 
489    {
490      strcpy(Value + strlen(Value), " ");
491    }
492
493  return elemidx;
494}
495
496
497int GET_GRIB_INDEX(FileIndex *fileindex, int *center, int *subcenter, 
498                   int *parmtbl, int *parmid, char DateStrIn[], 
499                   int *leveltype, int *level1, int *level2, int *fcsttime1,
500                   int *fcsttime2, int *index, int strlen1, int strlen2)
501{
502  char DateStr[1000];
503  FindGrib findgrib;
504
505  strncpy(DateStr,DateStrIn,strlen2);
506  DateStr[strlen2] = '\0';
507  grib_time_format(DateStr,DateStr);
508
509  rg_init_findgrib(&findgrib);
510
511  strncpy(findgrib.initdate,DateStrIn,strlen2);
512  findgrib.initdate[strlen2] = '\0';
513  findgrib.parmid          = *parmid;
514  findgrib.leveltype       = *leveltype;
515  findgrib.level1          = *level1;
516  findgrib.level2          = *level2;
517  findgrib.fcsttime1       = *fcsttime1;
518  findgrib.fcsttime2       = *fcsttime2;
519  findgrib.center_id       = *center;
520  findgrib.subcenter_id    = *subcenter;
521  findgrib.parmtbl_version = *parmtbl;
522 
523  *index = rg_get_index(fileindex->gribinfo, &findgrib);
524
525  return *index;
526
527}
528
529
530int GET_GRIB_INDEX_GUESS(FileIndex *fileindex, int *center, int *subcenter, 
531                         int *parmtbl, int *parmid, char DateStrIn[], 
532                         int *leveltype, int *level1, int *level2, 
533                         int *fcsttime1,int *fcsttime2, int *guessidx, 
534                         int *index, int strlen1, int strlen2)
535{
536  char DateStr[1000];
537  FindGrib findgrib;
538
539  strncpy(DateStr,DateStrIn,strlen2);
540  DateStr[strlen2] = '\0';
541  grib_time_format(DateStr,DateStr);
542
543  rg_init_findgrib(&findgrib);
544
545  strncpy(findgrib.initdate,DateStrIn,strlen2);
546  findgrib.initdate[strlen2] = '\0';
547  findgrib.parmid          = *parmid;
548  findgrib.leveltype       = *leveltype;
549  findgrib.level1          = *level1;
550  findgrib.level2          = *level2;
551  findgrib.fcsttime1       = *fcsttime1;
552  findgrib.fcsttime2       = *fcsttime2;
553  findgrib.center_id       = *center;
554  findgrib.subcenter_id    = *subcenter;
555  findgrib.parmtbl_version = *parmtbl;
556 
557  *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx);
558
559  return *index;
560
561}
562
563
564int GET_GRIB_CENTER(FileIndex *fileindex, int *parmid, int *center)
565{
566
567  *center = rg_get_center_id(fileindex->gribinfo,*parmid);
568
569  return *center;
570
571}
572
573
574int GET_GRIB_SUBCENTER(FileIndex *fileindex, int *parmid, int *subcenter)
575{
576
577  *subcenter = rg_get_subcenter_id(fileindex->gribinfo,*parmid);
578
579  return *subcenter;
580
581}
582
583
584int GET_GRIB_TBLVERSION(FileIndex *fileindex, int *parmid, int *parmtbl)
585{
586
587  *parmtbl = rg_get_parmtbl(fileindex->gribinfo,*parmid);
588
589  return *parmtbl;
590
591}
592
593
594int GET_GRIB_PROCID(FileIndex *fileindex, int *parmid, int *proc_id)
595{
596
597  *proc_id = rg_get_proc_id(fileindex->gribinfo,*parmid);
598
599  return *proc_id;
600
601}
602
603
604int GET_GRIB_INDEX_VALIDTIME(FileIndex *fileindex, int *center, 
605                   int *subcenter, int *parmtbl, int *parmid, 
606                   char DateStrIn[], int *leveltype, int *level1, int *level2,
607                   int *index, int strlen1, int strlen2)
608{
609  char DateStr[1000];
610  FindGrib findgrib;
611
612  strncpy(DateStr,DateStrIn,strlen2);
613  DateStr[strlen2] = '\0';
614  grib_time_format(DateStr,DateStr);
615
616  rg_init_findgrib(&findgrib);
617
618  strncpy(findgrib.validdate,DateStr,strlen2);
619  findgrib.initdate[strlen2] = '\0';
620  findgrib.parmid          = *parmid;
621  findgrib.leveltype       = *leveltype;
622  findgrib.level1          = *level1;
623  findgrib.level2          = *level2;
624  findgrib.center_id       = *center;
625  findgrib.subcenter_id    = *subcenter;
626  findgrib.parmtbl_version = *parmtbl;
627 
628  *index = rg_get_index(fileindex->gribinfo, &findgrib);
629
630  return *index;
631}
632
633
634int GET_GRIB_INDEX_VALIDTIME_GUESS(FileIndex *fileindex, int *center, 
635                                   int *subcenter, int *parmtbl, int *parmid, 
636                                   char DateStrIn[], int *leveltype, 
637                                   int *level1, int *level2, int *guessidx, 
638                                   int *index, int strlen1, int strlen2)
639{
640  char DateStr[1000];
641  FindGrib findgrib;
642
643  strncpy(DateStr,DateStrIn,strlen2);
644  DateStr[strlen2] = '\0';
645  grib_time_format(DateStr,DateStr);
646
647  rg_init_findgrib(&findgrib);
648
649  strncpy(findgrib.validdate,DateStr,strlen2);
650  findgrib.initdate[strlen2] = '\0';
651  findgrib.parmid          = *parmid;
652  findgrib.leveltype       = *leveltype;
653  findgrib.level1          = *level1;
654  findgrib.level2          = *level2;
655  findgrib.center_id       = *center;
656  findgrib.subcenter_id    = *subcenter;
657  findgrib.parmtbl_version = *parmtbl;
658 
659  *index = rg_get_index_guess(fileindex->gribinfo, &findgrib, *guessidx);
660
661  return *index;
662}
663
664
665int GET_GRIB_INDICES(FileIndex *fileindex, int *center, int *subcenter, 
666                     int *parmtbl,int *parmid, char DateStrIn[], 
667                     int *leveltype, int *level1, int *level2, int *fcsttime1,
668                     int *fcsttime2, int *indices, int *num_indices, 
669                     int strlen1, int strlen2)
670{
671  char DateStr[1000];
672  int status;
673  FindGrib findgrib;
674
675  strncpy(DateStr,DateStrIn,strlen2);
676  DateStr[strlen2] = '\0';
677  grib_time_format(DateStr,DateStr);
678
679  rg_init_findgrib(&findgrib);
680
681  strncpy(findgrib.initdate,DateStrIn,strlen2);
682  findgrib.initdate[strlen2] = '\0';
683  trim(findgrib.initdate);
684  findgrib.parmid          = *parmid;
685  findgrib.leveltype       = *leveltype;
686  findgrib.level1          = *level1;
687  findgrib.level2          = *level2;
688  findgrib.fcsttime1       = *fcsttime1;
689  findgrib.fcsttime2       = *fcsttime2;
690  findgrib.center_id       = *center;
691  findgrib.subcenter_id    = *subcenter;
692  findgrib.parmtbl_version = *parmtbl;
693
694  *num_indices = rg_get_indices(fileindex->gribinfo, &findgrib, indices);
695
696  return (*num_indices);
697
698}
699
700
701int GET_GRID_INFO_SIZE(int *size)
702{
703
704  *size = sizeof(Grid_Info);
705
706  return *size;
707
708}
709
710
711int LOAD_GRID_INFO(char *varnameIn, char *initdateIn, int *leveltype, 
712                   int *level1, int *level2, float *fcst_time, 
713                   int *accum_period, int *grid_id, int *projection, 
714                   int *xpoints, int *ypoints, float *center_lat, 
715                   float *center_lon, float *Di, float *Dj,float *central_lon,
716                   int *proj_center_flag, float *latin1, 
717                   float *latin2, Grib1_Tables *grib_tables,
718                   Grid_Info *grid_info, int strlen1, int strlen2)
719{
720
721  char varname[1000], initdate[1000];
722
723  strncpy(varname,varnameIn,strlen1);
724  varname[strlen1] = '\0';
725  strncpy(initdate,initdateIn,strlen2);
726  initdate[strlen2] = '\0';
727
728  strcpy(grid_info->varname, varname);
729  strcpy(grid_info->initdate, initdate);
730  grid_info->leveltype        = *leveltype;
731  grid_info->level1           = *level1 ;
732  grid_info->level2           = *level2 ;
733  grid_info->fcst_time        = *fcst_time ;
734  grid_info->accum_period     = *accum_period ;
735  grid_info->grid_id          = *grid_id ;
736  grid_info->projection       = *projection ;
737  grid_info->xpoints          = *xpoints ;
738  grid_info->ypoints          = *ypoints ;
739  grid_info->center_lat       = *center_lat ;
740  grid_info->center_lon       = *center_lon;
741  grid_info->Di               = *Di ;
742  grid_info->Dj               = *Dj ;
743  grid_info->central_lon      = *central_lon ;
744  grid_info->proj_center_flag = *proj_center_flag ;
745  grid_info->latin1           = *latin1 ;
746  grid_info->latin2           = *latin2 ;
747  grid_info->grib_tables      = copy_grib_tables(grib_tables);
748
749  return 0;
750
751}
752
753int PRINT_GRID_INFO(Grid_Info *grid_info)
754{
755
756  fprintf(stdout,"varname        =%s\n",grid_info->varname);
757  fprintf(stdout,"initdate       =%s\n",grid_info->initdate);
758  fprintf(stdout,"leveltype      =%d\n",grid_info->leveltype);
759  fprintf(stdout,"level1         =%d\n",grid_info->level1);
760  fprintf(stdout,"level2         =%d\n",grid_info->level2);
761  fprintf(stdout,"fcst_time      =%f\n",grid_info->fcst_time);
762  fprintf(stdout,"accum_period   =%d\n",grid_info->accum_period);
763  fprintf(stdout,"grid_id        =%d\n",grid_info->grid_id);
764  fprintf(stdout,"projection     =%d\n",grid_info->projection);
765  fprintf(stdout,"xpoints        =%d\n",grid_info->xpoints);
766  fprintf(stdout,"ypoints        =%d\n",grid_info->ypoints);
767  fprintf(stdout,"center_lat     =%f\n",grid_info->center_lat);
768  fprintf(stdout,"center_lon     =%f\n",grid_info->center_lon);
769  fprintf(stdout,"Di             =%f\n",grid_info->Di);
770  fprintf(stdout,"Dj             =%f\n",grid_info->Dj);
771  fprintf(stdout,"central_lon    =%f\n",grid_info->central_lon);
772  fprintf(stdout,"proj_center_flag =%d\n",grid_info->proj_center_flag);
773  fprintf(stdout,"latin1         =%f\n",grid_info->latin1);
774  fprintf(stdout,"latin2         =%f\n",grid_info->latin2);
775
776  return 0;
777
778}
779
780
781int GET_SIZEOF_GRID(FileIndex *fileindex, int *index, int *numcols, 
782                    int *numrows)
783{
784
785  *numcols = rg_get_numcols(fileindex->gribinfo,*index);
786
787  *numrows = rg_get_numrows(fileindex->gribinfo,*index);
788
789  return (*numcols)*(*numrows);
790
791}
792
793
794void FREE_GRID_INFO(Grid_Info *grid_info)
795{
796  FREE_GRIBMAP(grid_info->grib_tables);
797}
798
799
800int READ_GRIB(FileIndex *fileindex, int *fid, int *index, float *data)
801{
802  int status;
803
804  status = rg_get_data_1d(fileindex->gribinfo,*index,data);
805
806  return status;
807}
808
809#define LINESIZE 300
810
811#define SECS_IN_SEC 1
812#define SECS_IN_MIN 60
813#define MINS_IN_HOUR 60
814#define MINS_IN_5MINS 5
815#define HOURS_IN_DAY 24
816
817#define MAX_FCST 65535
818#define MAX_FCST_SECS MAX_FCST*SECS_IN_SEC
819#define MAX_FCST_MINS MAX_FCST*SECS_IN_MIN
820#define MAX_FCST_5MINS MAX_FCST*MINS_IN_5MINS*SECS_IN_MIN
821#define MAX_FCST_HOURS MAX_FCST*MINS_IN_HOUR*SECS_IN_MIN
822#define MAX_FCST_DAYS MAX_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN
823
824#define MAX1B_FCST 256
825#define MAX1B_FCST_SECS MAX1B_FCST*SECS_IN_SEC
826#define MAX1B_FCST_MINS MAX1B_FCST*SECS_IN_MIN
827#define MAX1B_FCST_5MINS MAX1B_FCST*MINS_IN_5MINS*SECS_IN_MIN
828#define MAX1B_FCST_HOURS MAX1B_FCST*MINS_IN_HOUR*SECS_IN_MIN
829#define MAX1B_FCST_DAYS MAX1B_FCST*HOURS_IN_DAY*MINS_IN_HOUR*SECS_IN_MIN
830
831typedef struct {
832  int time_range;
833  int fcst_unit;
834  int P1;
835  int P2;
836
837  int time_range_ext;
838  int fcst_unit_ext_1;
839  int fcst_unit_ext_2;
840  int P1_ext;
841  int P2_ext;
842} FcstTimeStruct;
843
844int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *fcst_time);
845
846/****************************************************************************
847 *
848 * This function takes in metadata in the grid_info structure, output data in
849 *   the *data array, and calls routines to write the metadata and data
850 *   in grib version 1 format the open file descriptor filefd.
851 *
852 ****************************************************************************/
853
854int WRITE_GRIB(Grid_Info *grid_info, int *filefd, float *data)
855{
856
857  GRIB_HDR *gh=NULL;
858  DATA_INPUT data_input;
859  GEOM_IN geom_in;
860  USER_INPUT user_input;
861  int grid_projection;
862  int status;
863  float x_center, y_center;
864  GridNav gridnav;
865  float first_lat, first_lon, last_lat, last_lon;
866  int year, month, day, hour, minute;
867  float second;
868  char varname2[1000];
869  int table_index;
870  int fcst_unit;
871  int time_range;
872  int P1, P2;
873  int fcst_unit_ext_1, fcst_unit_ext_2;
874  int P1_ext, P2_ext;
875  int time_range_ext;
876  char errmsg[1000];
877  int center, subcenter, parmtbl;
878  int tablenum;
879  float dx, dy;
880  FcstTimeStruct fcst_time;
881
882  strcpy(varname2,grid_info->varname);
883  trim(varname2);
884
885  sscanf(grid_info->initdate,"%d-%d-%d_%d:%d:%f",
886         &year,&month,&day,&hour,&minute,&second);
887
888  /* Get coords of center of grid */
889  x_center = (grid_info->xpoints + 1)/2.;
890  y_center = (grid_info->ypoints + 1)/2.;
891
892  grid_projection = get_gridnav_projection(grid_info->projection);
893
894 /* Convert grid spacing to degrees for LATLON projection */
895 
896  if ( (grid_info->projection == WRF_LATLON) || 
897       (grid_info->projection == WRF_CASSINI) )
898  {
899      dx = grid_info->Di * KM_TO_DEGREES;
900      dy = grid_info->Dj * KM_TO_DEGREES;
901  }
902  else
903  {
904      dx = grid_info->Di;
905      dy = grid_info->Dj;
906  }
907
908 /* Initialize grid structure */
909
910  status = GRID_init(grid_info->center_lat, grid_info->central_lon, 
911                     grid_projection,
912                     grid_info->latin1, grid_info->latin2, 
913                     grid_info->xpoints, grid_info->ypoints,
914                     dx, dy,
915                     grid_info->center_lat, grid_info->center_lon, 
916                     x_center, y_center,
917                     &gridnav);
918  if (!status)
919    {
920      fprintf(stderr,"write_grib: error from GRID_init\n");
921    }
922
923  /* get lat/lon of lower left corner */
924  status = GRID_to_latlon(&gridnav, 1, 1, &first_lat, &first_lon);
925  if (!status)
926    {
927      fprintf(stderr,
928              "write_grib: error from GRID_to_latlon for first lat/lon\n");
929    }
930
931  /* get lat/lon of upper right corner */
932  status = GRID_to_latlon(&gridnav, grid_info->xpoints, grid_info->ypoints, 
933                          &last_lat, &last_lon);
934  if (!status)
935    {
936      fprintf(stderr,
937              "write_grib: error from GRID_to_latlon for last lat/lon\n");
938    }
939
940  /* Read the grib parameter table */
941  status = GET_GRIB_PARAM(grid_info->grib_tables, varname2, &center,
942                           &subcenter, &parmtbl, &tablenum, &table_index,
943                          1,strlen(varname2));
944  if (table_index < 0)
945    {
946      fprintf(stderr,\
947              "Skipping %s, Could not find parameter for %s in gribmap.txt\n",\
948              varname2,varname2);
949      return 1;
950    }
951
952  /*
953   * We skip any parameters that are listed in parameter 255 in gribmap.txt.
954   * Parameter 255 is used to indicate that a WRF parameter should not be
955   * output.  It is useful for parameters that are requested to be output in
956   * the WRF Registry, but are already implicitly output in grib.
957   */
958
959  if (table_index == 255)
960    {
961      return 0;
962    }
963
964  /*
965   * Setup the geom_in structure for the grib library.  Here, we set
966   *  the generic parms.  Below, we set the projection specific parms
967   */
968  geom_in.nx = grid_info->xpoints;
969  geom_in.ny = grid_info->ypoints;
970  geom_in.first_lat = first_lat;
971  geom_in.first_lon = first_lon;
972  geom_in.last_lat = last_lat;
973  geom_in.last_lon = last_lon;
974  geom_in.scan = 64;
975
976  switch (grid_info->projection) 
977    {
978    case WRF_LATLON:
979    case WRF_CASSINI:
980      strcpy(geom_in.prjn_name,"spherical");
981      geom_in.parm_1 = dy;
982      geom_in.parm_2 = dx;
983      geom_in.parm_3 = -1;
984      geom_in.usRes_flag = 0;  /*
985                                * Set to 0 here, MEL grib library will reset
986                                *  to 128 to indicate that direction
987                                *  increments are given.
988                                */
989      break;
990    case WRF_MERCATOR:
991      strcpy(geom_in.prjn_name,"mercator");
992      geom_in.parm_1 = grid_info->latin1;
993      geom_in.parm_2 = dx;
994      geom_in.parm_3 = dy;
995      geom_in.usRes_flag = 136;
996      break;
997    case WRF_LAMBERT:
998      strcpy(geom_in.prjn_name,"lambert");
999      geom_in.usRes_flag = 0;   /* Set to 0 here, MEL grib library will reset
1000                                 * to 128.
1001                                 */
1002      geom_in.parm_3 = grid_info->central_lon;
1003      geom_in.x_int_dis = dx;
1004      geom_in.y_int_dis = dy;
1005      geom_in.parm_1 = grid_info->latin1;
1006      geom_in.parm_2 = grid_info->latin2;
1007      break;
1008    case WRF_POLAR_STEREO:
1009      strcpy(geom_in.prjn_name,"polar_stereo");
1010      geom_in.usRes_flag = 0;   /* Set to 0 here, MEL grib library will reset
1011                                 * to 128.
1012                                 */
1013
1014      geom_in.parm_3 = -1;
1015      geom_in.x_int_dis = dx*(1.+sin(60. * PI/180.))
1016        / (1.+sin(abs(grid_info->latin1) * PI/180.));
1017      geom_in.y_int_dis = dy*(1.+sin(60. * PI/180.))
1018        / (1.+sin(abs(grid_info->latin1) * PI/180.));
1019      geom_in.parm_1 = -1;
1020      geom_in.parm_2 = grid_info->central_lon;
1021      break;
1022    default:
1023      fprintf(stderr,"Error, invalid projection: %d\n",grid_info->projection);
1024      return 1;
1025    }
1026
1027  /*
1028   * Setup the data_input structure.
1029   */
1030  data_input.nDec_sc_fctr = 
1031    grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index];
1032  data_input.usProc_id = 220;
1033  data_input.usGrid_id = grid_info->grid_id;
1034  data_input.usParm_id = 
1035    grid_info->grib_tables->grib_table_info[tablenum].parm_id[table_index];
1036  data_input.usParm_sub_id = 
1037    grid_info->grib_tables->grib_table_info[tablenum].subcenter;
1038  data_input.usLevel_id = grid_info->leveltype;
1039
1040  if (grid_info->leveltype == 112) {
1041    data_input.nLvl_1 = grid_info->level2;
1042    data_input.nLvl_2 = grid_info->level1;
1043  } else {
1044    data_input.nLvl_1 = grid_info->level1;
1045    data_input.nLvl_2 = grid_info->level2;
1046  }
1047
1048  data_input.nYear = year;
1049  data_input.nMonth = month;
1050  data_input.nDay = day;
1051  data_input.nHour = hour;
1052  data_input.nMinute = minute;
1053  data_input.nSecond = second;
1054
1055  status = get_fcst_time(grid_info->accum_period, grid_info->fcst_time, 
1056                         &fcst_time);
1057
1058  data_input.usFcst_id = fcst_time.fcst_unit;
1059  data_input.usFcst_per1 = fcst_time.P1;
1060  data_input.usFcst_per2 = fcst_time.P2;
1061  data_input.usTime_range_id = fcst_time.time_range;
1062  data_input.usTime_range_avg = 0;
1063  data_input.usTime_range_mis = 0;
1064  /*
1065   * This is for WSI's extended PDS section
1066   */
1067  data_input.PDS_41 = fcst_time.fcst_unit_ext_1;
1068  data_input.PDS_42 = fcst_time.P1_ext;
1069  data_input.PDS_46 = fcst_time.fcst_unit_ext_2;
1070  data_input.PDS_47 = fcst_time.P2_ext;
1071  data_input.PDS_51 = fcst_time.time_range_ext;
1072  data_input.PDS_52 = 0;
1073
1074
1075  data_input.nDec_sc_fctr = 
1076    grid_info->grib_tables->grib_table_info[tablenum].dec_sc_factor[table_index];
1077  user_input.usCenter_id = 
1078    grid_info->grib_tables->grib_table_info[tablenum].center;
1079  user_input.usParm_tbl = 
1080    grid_info->grib_tables->grib_table_info[tablenum].parmtbl;
1081  user_input.chCase_id='0';
1082  user_input.usSub_tbl = 0;
1083  user_input.usCenter_sub = 
1084    grid_info->grib_tables->grib_table_info[tablenum].subcenter;
1085  user_input.usTrack_num = 0;
1086  user_input.usGds_bms_id = 128;
1087  user_input.usBDS_flag = 0;
1088  user_input.usBit_pack_num = 0;
1089
1090  status = init_gribhdr(&gh,errmsg);
1091  if (status != 0) {
1092    fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1093    return 1;
1094  }
1095
1096  status = grib_enc(data_input,user_input,geom_in,data,gh,errmsg);
1097  if (status != 0) {
1098    fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1099    fprintf (stderr,"\tCheck precision for %s in gribmap.txt.\n",
1100             varname2);
1101    return 1;
1102  }
1103
1104  status = gribhdr2filed(gh,*filefd,errmsg);
1105  if (status != 0) {
1106    fprintf (stderr,"write_grib: Error writing %s: \n\t%s\n",varname2,errmsg);
1107    return 1;
1108  }
1109
1110  free_gribhdr(&gh);
1111
1112  return 0;
1113}
1114
1115/***************************************************************************
1116 * Function to set up a structure containing forecast time parameters
1117 *  This encodes the standard grib forecast time parameters as well
1118 *  as WSI's extended forecast time parameters.
1119 ***************************************************************************/
1120
1121int get_fcst_time(int accum_period, int fcst_secs, FcstTimeStruct *ft)
1122{
1123  /*
1124   * Added ability to output a "5-minute" forecast time unit for the
1125   *   sake of WxProducer.  This allows WxPro to ingest data beyond
1126   *   18 hours, and accumulation data beyond 255 units.
1127   */
1128
1129  /*
1130   * Initialize.
1131   */
1132  ft->time_range = 0;
1133  ft->fcst_unit = 0;
1134  ft->P1 = 0;
1135  ft->P2 = 0;
1136  ft->fcst_unit_ext_1 = 0;
1137  ft->fcst_unit_ext_2 = 0;
1138  ft->P1_ext = 0;
1139  ft->P2_ext = 0;
1140  ft->time_range_ext = 0;
1141
1142  if (accum_period == 0) 
1143    {
1144      if (fcst_secs < MAX_FCST_SECS)
1145        {
1146          ft->time_range = 10;
1147          ft->fcst_unit = 254;
1148          ft->P1 = get_byte(fcst_secs,2);
1149          ft->P2 = get_byte(fcst_secs,1);
1150        }
1151      else if (((fcst_secs % SECS_IN_MIN) == 0) &&
1152               (fcst_secs < MAX_FCST_MINS))
1153        {
1154          ft->time_range = 10;
1155          ft->fcst_unit = 0;
1156          ft->P1 = get_byte(fcst_secs/SECS_IN_MIN,2);
1157          ft->P2 = get_byte(fcst_secs/SECS_IN_MIN,1);
1158        }
1159      else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR) == 0) && 
1160               (fcst_secs < MAX_FCST_HOURS))
1161        {
1162          ft->time_range = 10;
1163          ft->fcst_unit = 1;
1164          ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),2);
1165          ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR),1);     
1166        }
1167      /*
1168       * MAX_FCST_DAYS is causing an integer overflow, so, we'll just skip
1169       *   the check here.  It's very unlikely that someone would exceed this
1170       *   anyway (5.6 million days!)
1171       */
1172      /*
1173      else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0) &&
1174               (fcst_secs < MAX_FCST_DAYS))
1175      */
1176      else if (((fcst_secs % SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY) == 0))
1177        {
1178          ft->time_range = 10;
1179          ft->fcst_unit = 2;
1180          ft->P1 = 
1181            get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),2);
1182          ft->P2 = 
1183            get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR*HOURS_IN_DAY),1);
1184        }
1185      else if (((fcst_secs % SECS_IN_MIN*MINS_IN_5MINS) == 0)
1186               && (fcst_secs < MAX_FCST_5MINS))
1187        {
1188          ft->time_range = 10;
1189          ft->fcst_unit = 50;
1190          ft->P1 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),2);
1191          ft->P2 = get_byte(fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS),1);
1192        }
1193      else 
1194        {
1195          ft->time_range = 255;
1196          ft->fcst_unit = 0;
1197          ft->P1 = 0;
1198          ft->P2 = 0;
1199         
1200          ft->fcst_unit_ext_1 = 254;
1201          ft->fcst_unit_ext_2 = 254;
1202          ft->P1_ext = fcst_secs;
1203          ft->P2_ext = 0;
1204          ft->time_range_ext = 0;
1205        }
1206    }
1207  else  /* Accumulation period is not 0 */
1208    {
1209      if ((fcst_secs < MAX1B_FCST_HOURS) &&
1210          (fcst_secs%(SECS_IN_MIN*MINS_IN_HOUR) == 0) && 
1211          (accum_period%(SECS_IN_MIN*MINS_IN_HOUR) == 0))
1212        {
1213          ft->time_range = 4;
1214          ft->fcst_unit = 1;
1215          ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_HOUR);
1216          ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_HOUR);
1217        }
1218      else if ((fcst_secs < MAX1B_FCST_MINS) &&
1219               ((fcst_secs-accum_period)%SECS_IN_MIN == 0) && 
1220               (fcst_secs%SECS_IN_MIN == 0))
1221        {
1222          ft->time_range = 4;
1223          ft->fcst_unit = 0;
1224          ft->P1 = (fcst_secs-accum_period)/SECS_IN_MIN;
1225          ft->P2 = fcst_secs/SECS_IN_MIN;
1226        }
1227      else if (fcst_secs < MAX1B_FCST_SECS)
1228        {
1229          ft->time_range = 4;
1230          ft->fcst_unit = 254;
1231          ft->P1 = fcst_secs-accum_period;
1232          ft->P2 = fcst_secs;
1233        }
1234      else if ((fcst_secs < MAX1B_FCST_5MINS) && 
1235               (fcst_secs%(SECS_IN_MIN*MINS_IN_5MINS) == 0) &&
1236               (accum_period%(SECS_IN_MIN*MINS_IN_5MINS) == 0)) 
1237        {
1238          ft->time_range = 4;
1239          ft->fcst_unit = 50;
1240          ft->P1 = (fcst_secs-accum_period)/(SECS_IN_MIN*MINS_IN_5MINS);
1241          ft->P2 = fcst_secs/(SECS_IN_MIN*MINS_IN_5MINS);
1242        }
1243      else
1244        {
1245          ft->time_range = 255;
1246          ft->fcst_unit = 0;
1247          ft->P1 = 0;
1248          ft->P2 = 0;
1249         
1250          ft->fcst_unit_ext_1 = 254;
1251          ft->fcst_unit_ext_2 = 254;
1252          ft->P1_ext = fcst_secs - accum_period;
1253          ft->P2_ext = accum_period;
1254          ft->time_range_ext = 203;    /* Duration */
1255        }
1256    }
1257  return 0;
1258}
1259
1260
1261/******************************************************************************
1262 * returns a byt from an input integer
1263 *****************************************************************************/
1264
1265int get_byte(int input_int, int bytenum)
1266{
1267  int out;
1268  out = ((input_int >> (bytenum-1)*8) & ~(~0 <<8));
1269  return out;
1270}
1271
1272/*************************************************************************
1273 * Converts from WRF time format to time format required by grib routines
1274 *************************************************************************/
1275int grib_time_format(char *DateStr, char *DateStrIn)
1276{
1277  int year,month,day,hour,minute,second;
1278
1279  trim(DateStrIn);
1280  if (DateStrIn[0] == '*') {
1281    strcpy(DateStr,"*");
1282  }
1283  else
1284    {
1285      sscanf(DateStrIn,"%04d-%02d-%02d_%02d:%02d:%02d",
1286             &year,&month,&day,&hour,&minute,&second);
1287      sprintf(DateStr,"%04d%02d%02d%02d%02d%02d",
1288              year,month,day,hour,minute,second);
1289    }
1290
1291  return 0;
1292}
Note: See TracBrowser for help on using the repository browser.