ARJO SEGERS
KNMI
Date: 3 March 2009
Provides a module: use GribFile, only : ... The module provides a structure type containing all what should be known to open, read, and close a grib file: type(TGribFile) :: gribfile The gribfile is opened either for reading or writing with an 'Init' routine, which is in fact an interface to 'PbOpen'. A type TgribFile contains buffers to store one record of a grib file (a 'grib message'). 1. To open a grib file for reading, use: call Init( gribfile, filename, 'r', status ) type(TGribFile), intent(inout) :: gribfile character(len=*), intent(in) :: filename character(len=*), intent(in) :: mode integer, intent(out) :: status To read next record in buffer, use: call ReadRecord( gribfile, status ) Contents of the message is extracted with the 'Get' routine, for example: call Get( gribfile, status, pid=thepid ) Complete syntaxis and description of arguments (see also the 'GRIB SECTIONS' part later on): call Get( gribfile, status, & model_id, pid, & levtype, level, hyb_a, hyb_b, & reftime, timerange, & gridtype, & ll, & lon_first, lon_last, lon_inc, lon_n, & lat_first, lat_last, lat_inc, lat_n, & T, sh, & N, gg ) integer, intent(out), optional :: model_id ! Element 3 of section 1. integer, intent(out), optional :: pid ! Parameter identification. ! See the integer parameters 'pid_' in the code. ! Element 6 of section 1. integer, intent(out), optional :: levtype, level ! Elements 7-9 of section 1. ! --> Predefined values for level type: ! integer, parameter :: levtype_sfc = 1 ! surface ! integer, parameter :: levtype_hyb = 109 ! hybrid level real, intent(out), optional :: hyb_a(:), hyb_b(:) ! Vertical parameters in real section 2 integer, intent(out), optional :: reftime(5) ! Reference time: a 5 element integer array with ! year_inc_century, month, day, hour, minutes ! Elements 10-14,21 of section 1. integer, intent(out), optional :: timerange(4) ! Time range: a 4 element integer array with ! time range unit, value 1, value 2, indicator ! Elements 15-18 of section 1. integer, intent(out), optional :: gridtype ! Element 1 of section 2. ! --> Predefined values: ! integer, parameter :: gridtype_ll = 0 ! lat/lon ! integer, parameter :: gridtype_gg = 4 ! gaussian grid ! integer, parameter :: gridtype_sh = 50 ! spectral integer, intent(out), optional :: lon_first, lon_last, lon_inc, lon_n integer, intent(out), optional :: lat_first, lon_last, lat_inc, lat_n ! Elements of section 2 ! NOTE: values for lat/lon in mili degrees ! real, intent(out), optional :: ll(:,:) ! Value of lat/lon grid (stored from west->east, north->south!) complex, intent(out), optional :: sh(:) integer, intent(out), optional :: T ! spectral coefficients; ! size should match with spectral truncation 'T' real, intent(out), optional :: gg(:) integer, intent(out), optional :: N ! data values on Gaussian grid; ! number of values should match grid size 'N' Contents could be checked with the 'Check' routine, for example: call Check( gribfile, status, pid=154 ) Complete syntaxis: call Check( gribfile, status, debug=1, & model_id, pid, & levtype, level, & reftime, timerange, & gridtype, & lon_first, lon_last, lon_inc, lon_n, & lat_first, lat_last, lat_inc, lat_n, & T ) type(TGribFile), intent(in) :: gribfile integer, intent(out) :: status integer, intent(in), optional :: debug integer, intent(in), optional :: model_id integer, intent(in), optional :: pid integer, intent(in), optional :: levtype, level integer, intent(in), optional :: reftime(5) integer, intent(in), optional :: timerange(4) integer, intent(in), optional :: gridtype integer, intent(in), optional :: lon_first, lon_last, lon_inc, lon_n integer, intent(in), optional :: lat_first, lat_last, lat_inc, lat_n integer, intent(in), optional :: T integer, intent(in), optional :: status integer, intent(out), optional :: status See the 'Get' command above for a description. In case of one or more fields not matching the grib field: o return status is equal to number of tests failed; o info about the failed test is printed if debug is present and >0 . Note: hybride coefficients can not checked due to a lack of precission ... Similar for grid values (ll, gg, and sh). The grib file is closed with a 'Done' routine (interface to 'PbClose'): call Done( gribfile, status ) 2. To open a grib file for writing, use: call Init( gribfile, 'test.gb', 'w', status ) For safety, set all sections in the buffer to zero or an other safe value: call Clear( gribfile, status ) Fill some or all of the appropriate fields with the 'Set' command: call Set( gribfile, status, & model_id, pid, & levtype, hyb_a, hyb_b, level, & reftime, timerange, & gridtype, & ll, lon_n, lat_n, & lon_first, lon_inc, lon_last, & lat_first, lat_inc, lat_last, & scanning_mode, nbits ) type(TGribFile), intent(inout) :: gribfile integer, intent(out) :: status integer, intent(in), optional :: debug integer, intent(in), optional :: model_id integer, intent(in), optional :: pid integer, intent(in), optional :: levtype integer, intent(in), optional :: level real, intent(in), optional :: hyb_a(:), hyb_b(:) integer, intent(in), optional :: reftime(5) ! (/yy,mm,dd,hh,min/) integer, intent(in), optional :: timerange(4) integer, intent(in), optional :: gridtype real, intent(in), optional :: ll(:,:) integer, intent(in), optional :: lon_n, lat_n integer, intent(in), optional :: lon_first, lon_inc, lon_last ! mili degree integer, intent(in), optional :: lat_first, lat_inc, lat_last ! mili degree See the 'Get' command for a description of most of the fields. Only used as arguments to 'Set': integer, intent(in), optional :: scanning_mode ! Something with storage order; does not work proper yet. integer, intent(in), optional :: nbits ! Number of bits used to store a real value; default 24 . To write the record to the file, use: call WriteRecord( gribfile, status ) 3. To check wether a file is assigned to a grib structure, use the logical function 'Opened' : if ( Opened(gribfile) ) ... logical :: Opened type(TgribFile), intent(in) :: gribfile 4. To copy the buffer from one to another sturcture, use: call CopySections( gribIn, gribOut, status ) type(TGribFile), intent(in) :: gribIn type(TGribFile), intent(out) :: gribOut to copy all sections except real section 4 (the actual data); to copy the data too, use : call CopyAllSections( gribIn, gribOut, status ) type(TGribFile), intent(in) :: gribIn type(TGribFile), intent(out) :: gribOut
Compile together with the GribEx or Emos library from ECMWF:
f90 -o test.exe ... file_grib.o ... -l gribexSee also GRIBEX at the ECMWF website.
Overview of some useful entries in grib files. For code tables, see:
Section 1 - Product Definition Section. --------------------------------------- 1. Code Table 2 Version Number. 128 2. Originating centre identifier. 98 3. Model identification. 199 4. Grid definition. 255 5. Flag (Code Table 1) 10000000 6. Parameter identifier (Code Table 2). 7. Type of level (Code Table 3). 109 8. Value 1 of level (Code Table 3). 9. Value 2 of level (Code Table 3). 10. Year of reference time of data. 100 11. Month of reference time of data. 8 12. Day of reference time of data. 10 13. Hour of reference time of data. 12 14. Minute of reference time of data. 0 15. Time unit (Code Table 4). 1 16. Time range one. 0 17. Time range two. 0 18. Time range indicator (Code Table 5) 0 21. Century of reference time of data. Section 2 - Grid Description Section. ------------------------------------- Southern latitudes and Western longitudes are negative. Unit is mili degrees. 1. Data represent type = lat/long (Table 6) 0 2. Number of points along a parallel. 360 3. Number of points along a meridian. 181 4. Latitude of first grid point. 90000 5. Longitude of first grid point. -180000 6. Resolution and components flag. 10000000 7. Latitude of last grid point. -90000 8. Longitude of last grid point. 179000 9. i direction (East-West) increment. 1000 10. j direction (North-South) increment. 1000 11. Scanning mode flags (Code Table 8) 00000000 1. Data represent type = spectral (Table 6) 50 2. J - Pentagonal resolution parameter. 106 <--- T 3. K - Pentagonal resolution parameter. 106 4. M - Pentagonal resolution parameter. 106 5. Representation type (Table 9) 1 6. Representation mode (Table 10). 2 1. Data represent type = gaussian lat/lon (Table 6) 4 2. Number of points along a parallel <0=reduced, >0=regular 3. Number of points along a meridian. 160 4. Latitude of first grid point. 89141 5. Longitude of first grid point. 0 6. Resolution and components flag. 00000000 7. Latitude of last grid point. -89141 8. Longitude of last grid point. 358878 9. i direction (East-West) increment Not given 10. Number of parallels between pole and equator. 80 <--- N 11. Scanning mode flags (Code Table 8) 00000000 23-22*2N. Number of points along a parallel for reduced grid 12. Number of vertical coordinate parameters. 122 (parameters stored in real part of section 2, starting form index 11) Section 4 - Binary Data Section. ------------------------------------- 1. Number of data values coded/decoded. 11556 2. Number of bits per data value. 16 3. Type of data (0=grid pt, 128=spectral). 128 4. Type of packing (0=simple, 64=complex). 64 5. Type of data (0=float, 32=integer). 0 6. Additional flags (0=none, 16=present). 0 7. Reserved. 0 8. Number of values (0=single, 64=matrix). 0 9. Secondary bit-maps (0=none, 32=present). 0 10. Values width (0=constant, 16=variable). 0 11. Byte offset of start of packed data (N). 2214 12. Power (P * 1000). 500 13. Pentagonal resolution parameter J for subset. 20 14. Pentagonal resolution parameter K for subset. 20 15. Pentagonal resolution parameter M for subset. 20 3. Type of data (0=grid pt, 128=spectral). 0 4. Type of packing (0=simple, 64=complex). 0 5. Type of data (0=float, 32=integer). 0 6. Additional flags (0=none, 16=present). 0 7. Reserved. 0 8. Number of values (0=single, 64=matrix). 0 9. Secondary bit-maps (0=none, 32=present). 0 10. Values width (0=constant, 16=variable). 0