ecCodes 学习 利用ecCodes fortran90 api对GRIB文件进行读写

参考 https://www.ecmwf.int/assets/elearning/eccodes/eccodes2/story_html5.html
https://confluence.ecmwf.int/display/OPTR/ecCodes%3A+GRIB+data+decoding+and+encoding+software+2018

基本解码流程

1. 指定打开方式（“读”或“写”），打开一个或多个GRIB文件;

2. 根据不同加载方式，加载一个或多个GRIB messages到内存：
　　有两种函数：codes_grib_new_from_file 和 codes_new_from_index。调用后会返回一个唯一的identifier，用于对已加载的GRIB messages进行操纵。

3. 调用codes_get函数对已加载的GRIB messages进行解码; (可以解码需要的数据)

4. 释放已经加载的GRIB messages：
　　codes_release　　

5. 关闭打开的 GRIB 文件.

此外，eccodes还有以下功能：

　　eccodes的主要目的是提供一个高水平的方法，用以从一个加载的GRIB messages对象去提取和计算出更多额外的信息：

· 返回平均，最小，最大，特定的经纬度等的关键字； 

· 计算经纬度和值：codes_grib_get_data； 

· 提取值的子例程：
　　codes_grib_find_nearest 去提取出距离给定地理点最近位置的点的值
　　codes_get_element 从列表中提取值

· 基于索引访问(indexed access)的子例程：这些是随机访问的方法，通常比顺序访问更快。

 

-----------------------------------------------------------------

顺序访问方式：

大致思路：

-> codes_open_file 

-> codes_grib_new_from_file -> codes_get -> codes_release

     …

-> codes_grib_new_from_file -> codes_get-> codes_release 

-> codes_close_file

 

示例代码：

 

 

索引访问方式（通常比顺序访问快）：

注意，eccodes中的index文件(后缀为.idx)与GrADS中后缀为.idx的文件不能通用！

大致思路：

-> codes_index_create(从grib文件创建index) 或 codes_index_read（读取已有index）

-> codes_index_select 选取键值

-> codes_new_from_index -> codes_get -> codes_release

     …

-> codes_new_from_index -> codes_get -> codes_release

-> codes_index_release

 

示例代码：

-----------------------------------------------------------------

一些ecCodes命令行工具：

参考 https://confluence.ecmwf.int/display/GRIB/GRIB+tools

grib_filter 筛选

grib_filter [options] rules_file grib_file gribfile ...

options

-f　　强制执行

-o　　输出index文件。如果没指定输出文件，那么输出GRIB文件写在filtered.out

-M　　关闭多场支持。关闭在单GRIB message中多要素场的支持。

-V　　版本

-g　　复制GTS头

-G　　兼容GRIBEX模式

-T T|B　　message类型。T -> GTS, B -> BUFR, 输入文件类型根据messgae解释

-7　　当message长度错误时候不报错

-v　　verbose模式

例子

1. grib_filter程序对于输入文件中的所有GRIB message顺序处理，对其中各项均应用rules。输入消息可以使用“write”声明写入输出文件。write语句可以parameterised以便输出发送到多个文件，根据键值定义输出文件的名字。如果我们写一个rules_file包含唯一的声明：

write "../data/split/[centre]_[dataDate]_[dataType]_[levelType].grib[edition]";

将这个rules_file应用到 "../data/tigge_pf_ecmwf.grib2"文件将会在 ../data/split 目录下得到几个文件，包含根据键值分割的字段

> grib_filter rules_file ../data/tigge_pf_ecmwf.grib2
> ls ../data/split
ecmf_20060619_pf_sfc.grib2
ecmf_20060630_pf_sfc.grib2
ecmf_20070122_pf_pl.grib2
ecmf_20070122_pf_pt.grib2
ecmf_20070122_pf_pv.grib2
ecmf_20070122_pf_sfc.grib2

2. 通过明确指示冒号后所需的类型，也可以以不同的格式获取文件名中的键值。

• :s　　整型
• :d　　双精度
• :s　　字符串

以下语句的工作方式与前一个示例略有不同，包括输出文件名中的center和dataType的整数值。

write "../data/split/[centre:i]_[dataDate]_[dataType:i]_[levelType].grib[edition]";

再次运行相同的命令，我们获得了不同的文件列表。

> grib_filter rules_file ../data/tigge_pf_ecmwf.grib2
> ls ../data/split
98_20060619_4_sfc.grib2
98_20060630_4_sfc.grib2
98_20070122_4_pl.grib2
98_20070122_4_pt.grib2
98_20070122_4_pv.grib2
98_20070122_4_sfc.grib2

3. grib_filter语法中允许使用其他语句：

• 以＃开头的注释
• print "string to print also with key values like in the file name"
  
• transient keyname1 = keyname2;
• set keyname = keyvalue;
• 定义（keyname）以检查message中是否定义了关键字
• 缺少（keyname）来检查关键字的值是否设置为MISSING
• 要将键值设置为MISSING，请使用“set key = MISSING;” （见示例）
• if（condition）{block of rules} else {block of rules} 
  条件可以使用==，！=，并使用||和&& 连接单个块条件 
  语句可以是任何有效语句也是另一个嵌套条件
• 您也可以使用“assert（condition）”进行断言。如果condition为false，它将中止filter。
  例如：assert（edition == 1）;

grib_filter规则的一个复杂示例如下是在GRIB版本1文件中更改温度。

# This filter should only be run on GRIB edition 1; abort otherwise
assert( edition == 1 );
 
# Temperature
if ( level == 850 && indicatorOfParameter == 11 ) {
    print "found indicatorOfParameter=[indicatorOfParameter] level=[level] date=[dataDate]";
    transient oldtype = type ;
    set identificationOfOriginatingGeneratingSubCentre=98;
    set gribTablesVersionNo = 128;
    set indicatorOfParameter = 130;
    set localDefinitionNumber=1;
    set marsClass="od";
    set marsStream="kwbc";
    # Negatively/Positively Perturbed Forecast
    if ( oldtype == 2 || oldtype == 3 ) {
      set marsType="pf";
      set experimentVersionNumber="4001";
    }
    # Control Forecast
    if ( oldtype == 1 ) {
      set marsType="cf";
      set experimentVersionNumber="0001";
    }
    set numberOfForecastsInEnsemble=11;
    write;
    print "indicatorOfParameter=[indicatorOfParameter] level=[level] date=[dataDate]";
    print;
}

4. 以下是将键与字符串进行比较的IF语句示例。注意，必须对字符串使用“is”关键字而不是“==”。并且，为了否定，您需要在整个条件之前添加“！”：

# Select Geopotential Height messages which are not on a Reduced Gaussian Grid
if (shortName is "gh" && !(gridType is "reduced_gg" )) {
    set step = 72;
}

 

5. switch语句是if语句的替代版本。其语法如下：

switch (key1) {
    case val1:
        # block of rules;
    case val2:
        # block of rules;
    default:
        # [ block of rules ]
}

 

作为switch语句的参数给出的键的每个值都与case语句中指定的值匹配。如果存在匹配，则执行与匹配的case语句对应的块或规则。否则，执行默认情况。如果case声明没有包含所有可能性，则默认情况是强制性的。“〜”运算符可用于匹配“任何”。

以下示例显示了switch语句的用法：

processing paramId=[paramId] [shortName] [stepType]
switch (shortName) {
    case "tp":
        set stepType="accum";
    case "10u" :
        set typeOfLevel="surface";
  default:
}

  

grib_get 获取数据

Get values of some keys from a GRIB file. It is similar to grib_ls, but fails returning an error code when an error occurs (e.g. key not found).

grib_get [options] grib_file grib_file ...

options

-f　　强制执行

-p key[:{s/d/i}],key[:{s/d/i}],...
　　声明打印的关键字。可以要求字符串型(keys:s)，双精度(keys:d)或者整型(keys:i)的关键字。默认类型是字符串

-F format　　C样式格式的浮点值。

-l Latitude,Longitude[,MODE,file]
　　　接近的纬度/经度点的值。允许的MODE有：

• 4 　　默认 (打印距离点最近四个点的值)
• 1 　　打印距离点最近的值
• file 　文件被用于mask。打印距离mask值>=0.5时候，最近的点的值

-P key[:{s/d/i}],key[:{s/d/i}],...　　As -p adding the declared keys to the default list.

-w key[:{s/d/i}]{=/!=}value,key[:{s/d/i }]{=/!=}value,...
　　只有满足所有设置键/值约束的GRIB message才会被处理。一个有效的键值约束类似于key=value或者key!=value。对于任意一个关键字，可以是字符串型(keys:s)，双精度型(keys:d)或者一个长整型(keys:l)。默认类型为字符串型。在值中同样可以使用上斜杠“/"去指定一个 或 条件。(即一个逻辑析取)

-n namespace　　打印所有属于namespace的关键字。一些有用的namespace是:“time”,“parameter”,“geography”和“statistics”。

-s key[:{s/d/i}]=value,key[:{s/d/i}]=value,...
　　　　需要设置的关键字数值。对于每一个关键字，可以定义为字符串(key:s)、双精度(key:d)或者整型(key:i)这些类型。默认情况是设置为native type

-V　　Version. 

-W width　　输出行的最小宽度。默认是10

-M　　关闭多场支持。关闭在单GRIB message中多要素场的支持。

-g　　复制GTS头

-G　　兼容GRIBEX模式

-T T|B　　message类型。T -> GTS, B -> BUFR, 输入文件类型根据messgae解释

-7　　当message长度错误时候不报错

-i index　　打印对应于给定index的值。注意这个index是从零开始的，所以第一个值是在索引0。

例子

1. 如果关键字没有找到，grib_get 失败

> grib_get -p gribname ../data/tigge_pf_ecmwf.grib2

2. 为了获取文件中第一个GRIB message中step的值 

> grib_get -w count=1 -p step ../data/tigge_pf_ecmwf.grib2

  

grib_index_build 创建索引 

grib_index_build [options] grib_file grib_file ...

options 

-f　　强制执行

-o　　输出index文件。如果没指定输出文件，那么输出文件写成“gribidx”

-k　　keys1,keys2... 给定需要索引的关键字。默认输入文件将用MARS关键字索引。可以要求字符串型(keys:s)，双精度(keys:d)或者整型(keys:i)的关键字。

-V　　版本号

-M　　多变量场支持关闭。关闭单独GRIB message中多变量场的支持

-N　　不要压缩索引。默认情况下，索引移除了只有一个值的关键字，用以压缩。

例子

1. 默认情况下，grib_index将会用MARS关键字进行索引

> grib_index_build ../data/reduced*.grib1 ../data/regular*.grib1 ../data/reduced*.grib2

2. 如果要给索引设置默认的关键字，使用-k选项。

> grib_index_build ../data/reduced*.grib1 ../data/regular*.grib1 ../data/reduced*.grib2

注意！因为并非所有keys都会被保存在grib文件中，一些keys是通过其它keys计算得到的。而在利用grib_index_build命令时候，默认只存储部分keys，如果需要在index中保存额外的keys，还需要打开-N选项：

$ grib_index_build -k time,name,level,longitude,latitude,month,day,year,shortName -N  *.grb -o outindex
--- grib_index_build: processing t.20180831.grb
--- grib_index_build: keys included in the index file outindex:
--- time, name, level, longitude, latitude, month, day, year, shortName
--- time = { 0, 600, 1200, 1800 }
--- name = { Temperature }
--- level = { 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70 }
--- longitude = { undef }
--- latitude = { undef }
--- month = { 8 }
--- day = { 31 }
--- year = { 2018 }
--- shortName = { t }
--- 112 messages indexed

如果不加-N选项，那么一些用户定义的keys就没有保存到index文件中

$ grib_index_build -k time,name,level,longitude,latitude,month,day,year,shortName   *.grb -o outindex
--- grib_index_build: processing t.20180831.grb
--- grib_index_build: keys included in the index file outindex:
--- time, level
--- time = { 0, 600, 1200, 1800 }
--- level = { 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70 }
--- 112 messages indexed

如上图，一些变量就没有保存到index文件中。

 

grib_dump 显示一个索引文件的内容

grib_dump [options] grib_file grib_file ...

options 

-O　　Octet（八进制？）模式。WMO文档形式的输出

-D　　调试模式

-d　　仅在C模式下可用。打印所有数据值。

-C　　C编码模式。一个C编码程序生成的GRIB message被输出

-t　　打印类型信息

-H　　以十六位格式打印八位内容

-a　　打印别名（aliases）

-w key[:{s/d/l}]{=/!=}value,key[:{s/d/l}]{=/!=}value,...  
　　　　只有满足所有设置键/值约束的GRIB message才会被处理。一个有效的键值约束类似于key=value或者key!=value。对于任意一个关键字，可以是字符串型(keys:s)，双精度型(keys:d)或者一个长整型(keys:l)。默认类型为字符串型。

-s key[:{s/d/i}]=value,key[:{s/d/i}]=value,...
　　　　需要设置的关键字数值。对于每一个关键字，可以定义为字符串(key:s)、双精度(key:d)或者整型(key:i)这些类型。默认情况是设置为native type

-M　　关闭多场支持。关闭在单GRIB message中多要素场的支持。

-T T|B|A　　message类型。T -> GTS, B -> BUFR, A -> Any (Experimental). 输入文件类型根据messgae解释

-7　　当message长度错误时候不报错

-V　　版本

-X offset　　输入文件offset量（单位字节）。处理输入文件将从“offset”开始。

例子

1. 用hexadecimal octet的方式(-H)，用WMO文档的形式输出。

> grib_dump -OH ../data/reduced_gaussian_model_level.grib1

2. 添加关键字别名和类型信息

> grib_dump -OtaH ../data/reduced_gaussian_model_level.grib

3. 获得一个grib文件中所有的关键字名称(包括已经计算得到的关键字)

> grib_dump -D ../data/regular_latlon_surface.grib1

  

grib_copy 复制GRIB文件的内容打印一些键的值 

grib_copy [options] grib_file grib_file ... output_grib_file

options 

-f　　强制执行

-r 　　重新包装数据。有时在设置一些涉及打包算法属性的键之后，需要重新打包数据。执行此重新打包时，设置此-r选项

-p key[:{s/d/l}],key[:{s/d/l}],... 
Declaration of keys to print. For each key a string (key:s) or a double (key:d) or a long (key:l) type can be requested. Default type is string. 

-P key[:{s/d/l}],key[:{s/d/l}],... 
As -p adding the declared keys to the default list. 

-w key[:{s/d/l}]{=/!=}value,key[:{s/d/l}]=value,... 
　　　　只有满足所有设置键/值约束的GRIB message才会被处理。一个有效的键值约束类似于key=value或者key!=value。对于任意一个关键字，可以是字符串型(keys:s)，双精度型(keys:d)或者一个长整型(keys:l)。默认类型为字符串型。
　　　　Where clause. Only grib messages matching the key/value constraints are copied to the output_grib_file.
　　　　A valid constraint is of type key=value or key!=value.
　　　　For each key a string (key:s) or a double (key:d) or a long (key:l) type can be defined. Default type is string.
　　　　In the value you can also use the forward-slash character "/" to specify an OR condition (i.e. a logical disjunction)

-B order by directive 
Order by. The output will be ordered according to the "order by" directive. For example: "step:i asc, centre desc" (step numeric ascending and centre descending) 
-V　　版本号
-W width 
Minimum width of each column in output. Default is 10. 

-M 　　多变量场支持关闭。关闭单独GRIB message中多变量场的支持 Multi-field support off. Turn off support for multiple fields in single grib message.

-T T | B
Message type. T->GTS, B->BUFR. The input file is interpreted according to the message type.

-g 
Copy GTS header. 

-G 
GRIBEX compatibility mode. 

-7 
Does not fail when the message has wrong length 

-v 
Verbose.

 

例子

1. To copy only the pressure levels from a file

> grib_copy -w levtype=pl ../data/tigge_pf_ecmwf.grib2 out.grib

 

2.To copy only the fields that are not on pressure levels from a file

> grib_copy -w levtype!=pl ../data/tigge_pf_ecmwf.grib2 out.grib

3.To copy only the first three fields from a file

> grib_copy -w count=1/2/3 ../data/tigge_pf_ecmwf.grib2 out.grib

 

4. A grib_file with multi field messages can be converted in single field messages with a simple grib_copy

> grib_copy multi.grib simple.grib

 

5. Use the square brackets to insert the value of a key in the name of the output file (This is a good way to split a large GRIB file)
Note: we need to quote the name of the output so the shell does not interpret the square brackets

> grib_copy in.grib 'out_[shortName].grib'

 

6. To copy fields whose typeOfLevel is either "surface" or "meanSea"

> grib_copy -w typeOfLevel=surface/meanSea orig.grib out.grib

 

7. To copy selected fields and apply sorting (sorted by level in ascending order)
Note: we need to specify the ":i" to get a numerical sort. By default values are sorted as strings so a level of 100 would come before 20!

> grib_copy -w typeOfLevel=heightAboveGround -B"level:i asc" tigge_af_ecmwf.grib2 out.grib

-----------------------------------------------------------------

部分函数参考：

codes_get (msgid, key, value, status)
　　Get the value for a key from a grib message
　　从grib message中获取键值。
　　输入一个msgid和key，返回关键字的值。 在某些情况下，值可以是数组而不是标量。 作为数组键的示例，我们分别具有“values”，“pl”，“pv”数据值，简化网格中每个纬度的点数列表以及垂直级别列表。 在这些情况下，值数组必须由调用者分配，并且可以使用codes_get_size获取其所需的维度。值可以是整数（4），实数（4），实数（8），字符。 虽然每个键都有自己的本机类型，但是可以检索整数类型的键（使用codes_get）作为real（4），real（8）或character。 尽可能提供类似的转换。 对于任何其他类型，非法转换对整数和字符都是真实的。msgid引用内存中加载的消息。
　　[in] msgid id of the message loaded in memory
　　[in] key key name
　　[out] value value can be a scalar or array of integer(4),real(4),real(8),character. Arrays must support the allocatable attribute.
　　[out] status CODES_SUCCESS if OK, integer value on error

codes_get_element (msgid, key, kindex, value, status)
　　Get a value of specified index from an array key
　　从一个关键字数组中获取给定索引的值。给定 message id、关键字名称、和索引，得到相应的值。 索引是从0开始的（即第一个元素的索引是0，第二个元素索引为1，依此类推）。 如果索引参数是数组，则返回与索引数组对应的所有值。
　　[in] id　　 加载到内存的ID。 ID of the message loaded in memory
　　[in] key 　关键字名称　key name
　　[in] index 　可以是一个标量或者数组。单精度整型　index can be a scalar or array of integer(4)
　　[out] value 　可以是一个标量或者数组。单精度整型，单精度或双精度实型数组。　value can be a scalar or array of integer(4),real(4),real(8)
　　[out] status 　　CODES_SUCCESS if OK, integer value on error

codes_get_error_string (error, error_message, status)
　　Get the error message given an error code.
　　输入错误码，返回错误消息
　　error　　　　　　　错误码
　　error messages 　　错误消息
　　status　　　　　　如果成功返回CODES_SUCCESS，错误返回整数值

codes_get_message_size (msgid, nbytes, status)
　　Get the size of a coded message
　　得到编码message的大小

codes_get_size (msgid, key, size, status)
　　Get the size of an array key
　　得到键数组的大小

codes_grib_find_nearest_single (gribid, is_lsm, inlat, inlon, outlat, outlon, value, distance, kindex, status)

　　获得距给定经纬度点最近的点的位置及值。
　　输入gribid，逻辑值
　　返回最近的一个点 (或最近的四个点) 的值，基于零的索引 (可在 code_get_element中使用) ，及其与给定点的距离。使用以下公式： radius * acos( sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2) ).
　　如果is_lsm flag是. true。输入场网格被认为是land sea mask, 返回最近的陆地点。
　　在四个相邻点中, 最近的陆地点是:
　　　　land sea mask值 > = 0.5 的最近点
　　　　如果这四个相邻点的land sea mask 值都 < 0.5, 则在没有任何其他条件的情况下最近的。
　　可以提供经纬度的数组 (real(8), 以便通过一个调用来查找数组中列出的所有 lat/lon 点的值、索引和距离。
　　如果提供了单个纬度/经度点, 并且outlat, outlon, value, 距离, 索引被定义为具有四个元素的数组, 返回四个最近点的 lat/lon 坐标和value、距离和索引。
　　如果出现错误, 如果未给出状态参数 (可选), 程序将退出并显示错误消息。否则, 可以使用codes_get_error_string 收集错误消息。
　　[in] gribid　　加载到内存中的值
　　[in] is_lsm　　逻辑值，如果要求最近的点，则为.true.，否则为.false.
　　[in] inlat　　点的纬度
　　[in] inlon　　点的经度
　　[out] outlat　　最近点的纬度
　　[out] outlon　　最近点的经度
　　[out] distance　　给点和最近点的距离
　　[out] kindex　　基于0的索引
　　[out] value　　最近点的值
　　[out] status　　状态 CODES_SUCCESS if OK, integer value on error 

codes_grib_get_data (gribid, lats, lons, values, status)
　　Get latitude/longitude and data values. 
　　得到经纬度和数据值

codes_index_add_file (indexid, filename, status)
　　Add a file to an index.
　　将一个文件加入到索引文件中

codes_index_create (indexid, filename, keys, status)
　　to create the index of the content of a file 创建文件内容索引
　　indexid　　新建索引文件的id
　　filename　　被索引messages的文件的名称
　　keys　　　　用逗号分隔的关键字关键字的类型可以在后面追加(:l，表示长整型，:i短整型，:d双精度，:s字符串）如果没有显式声明类型，那么假设类型默认。
　　status　　　如果成功返回CODES_SUCCESS，错误返回整数值

codes_index_get (indexid, key, values, status)
　　Get the distinct values of the key in argument contained in the index. 
　　得到包含在索引中声明的特定关键字的特定值

codes_index_get_size (indexid, key, size, status)
　　to get the dimension of a key in the index
　　得到索引的关键字的维数

codes_index_read (indexid, filename, status)
　　Load an index file previously created with codes_index_write.
　　加载一个索引文件（该索引文件之前已由codes_index_write命令创建）

codes_index_release (indexid, status)
　　释放index
　　indexid 　　id of an index created from a file.
　　status 　　 CODES_SUCCESS if OK, integer value on erro

codes_index_select (indexid, key, value, status)
　　Select the message subset with key==value.
　　用key==valued选取消息子集

codes_index_write (indexid, filename, status)
　　Saves an index to a file for later reuse.
　　保存索引文件便于以后使用

codes_new_from_file (ifile, msgid, product_kind, status)

　　从该文件中加载一个message 到内存。输入文件id、数据类型，并返回message id
　　可以通过msgid直接访问该message，msgid一直可用，调用codes_release函数之前。 

codes_new_from_index (indexid, msgid, status)
　　在选定key值后，从一个index创建一个新的handle。
　　在调用这个函数前，index的所有keys都比必须确定。 连续调用这个函数 将返回index keys中 所有符合定义的handles。当没有更多的可用处理从索引中返回一个空指针变量，并且err变量值设置为CODES_END_OF_INDEX。

codes_new_from_message(msgid, message, status)
　　Create a new message in memory from an integer or character array containting the coded message. 在包含编码消息的整型或字符型数组 中在内存中创建新消息

codes_grib_new_from_samples (gribid, samplename, status)
　　Create a new valid gribid from a GRIB sample contained in a samples directory pointed by the environment variable ECCODES_SAMPLES_PATH.
　　在ECCODES_SAMPLES_PATH环境变量指定的样本路径中的GRIB样本，创建新的可用的gribid.(编码时候用)
　　gribid 　　id of the grib loaded in memory 输出量，加载到内存中的message的id
　　samplename　　 name of the sample to be used 输入量，样本名称
　　status 　　CODES_SUCCESS if OK, integer value on error 输出量，是否成功

codes_open_file (ifile, filename, mode, status)
　　Open a file according to a mode. 根据模式打开文件
　　ifile 　　id of the opened file to be used in all the file functions.
　　filename 　　name of the file to be open
　　mode 　　open mode can be 'r' (read only), 'w' (write only) or 'a' (append)
　　status 　　如果成功返回CODES_SUCCESS，错误返回整数值

codes_read_bytes (ifile, buffer, nbytes, status)
　　Reads nbytes bytes into the buffer from a file opened with codes_open_file

codes_read_from_file (ifile, buffer, nbytes, status)
　　Reads a message in the buffer array from the file opened with codes_open_file.

codes_release (msgid, status)
　　Free the memory for the message referred as msgid.

codes_set (msgid, key, value, status)
　　Set the value for a key in a message.
　　在message中设置键值

codes_write_bytes_int4 (ifile, buffer, nbytes, status)
　　Write nbytes bytes from the buffer in a file opened with codes_open_file.

 注 ：ERROR CODES列表参考 http://download.ecmwf.int/test-data/eccodes/html/group__errors.html

 

编译示例

引用时候用-leccodes_f90这个库

假设存放eccodes的头文件和库文件路径的环境变量分别是ECCODES_INCLUDE和ECCODES_LIB

jlz@dell:~/test$ ifort test.f90 -I$ECCODES_INCLUDE -L$ECCODES_LIB -leccodes_f90 -o test.out

 

 

#2020-11-27更新#

用ECcodes写GRIB文件

大致流程

codes_grib

 

参考https://confluence.ecmwf.int/display/ECC/grib_set_data

! Copyright 2005-2018 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!
!  Description: set the data contained in a GRIB file.
!               In this example no missing values are present.
!               If there are missing values, refer to: grib_set_bitmap
!
program set_data
  use eccodes
  implicit none
  integer                         :: outfile !输出文件的id
  integer                         :: i, igrib, iret, numberOfValues, cnt
  real                            :: d, e
  real, dimension(:), allocatable :: values
  integer, parameter              :: max_strsize = 200
  character(len=max_strsize)      :: outfile_name
 
  call getarg(1, outfile_name) !执行的时候从外界获取参数作为outfile_name
  call codes_open_file(outfile,outfile_name,'w')  !以写模式打开输出文件，其id 为outfile
 
  ! Note: the full name of the sample file is "regular_ll_pl_grib1.tmpl"
  ! Sample files are stored in the samples directory (use codes_info to
  ! see where that is).  The default sample path can be changed by
  ! setting the environment variable ECCODES_SAMPLES_PATH
  call codes_grib_new_from_samples(igrib, 'regular_ll_pl_grib1')
 
  ! Here we're changing the data values only, so the number of values
  ! will be the same as the sample GRIB.
  ! But if your data array has a different size, then specify the grid geometry
  ! (e.g. keys Ni, Nj etc) and set the correct number of data values
  call codes_get_size(igrib,'values',numberOfValues)
 
  allocate(values(numberOfValues), stat=iret)
  d = 10e-8
  e = d
  cnt = 1
  do i=1,numberOfValues
    if (cnt>100) then
      e = e*10
      cnt=1
    endif
    values(i) = d
    !print *, values(i)
    d = d + e
    cnt = cnt + 1
  end do
 
  call codes_set(igrib, 'bitsPerValue', 16)
 
  ! set data values
  call codes_set(igrib,'values', values)
  call codes_write(igrib,outfile)
  call codes_release(igrib)
  deallocate(values)
 
end program set_data

 

grib_samples https://confluence.ecmwf.int/display/ECC/grib_samples

! Copyright 2005-2018 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!  Description: how to create a new GRIB message from a sample.
!   从例子中创建新的GRIB消息
!
program sample
  use eccodes
  implicit none
  integer  :: err
  integer  :: outfile, datafile
  integer  :: igribsample,igribclone,igribdata, size1
  integer  :: date1, startStep, endStep, table2Version, indicatorOfParameter
  integer  :: decimalPrecision
  character(len=10) stepType
  real(kind=8), dimension(:), allocatable   :: v1,v2,v
 
  date1 = 20080104
  startStep = 0
  endStep = 12
  stepType = 'accum'
  table2Version = 2
  indicatorOfParameter = 61
  decimalPrecision = 2
 
  ! A new grib message is loaded from an existing sample.
  ! Samples are searched in a default sample path (use codes_info
  ! to see where that is). The default sample path can be changed by
  ! setting the environment variable ECCODES_SAMPLES_PATH

  !从示例文件(存放在 ECCODES_samples_PATH指向的目录中)创建新

  call codes_grib_new_from_samples(igribsample, "regular_latlon_surface.grib1") 
  
  call codes_open_file(outfile, 'f_out.samples.grib1','w')   ! 创建输出文件

  call codes_open_file(datafile,'../../data/tp_ecmwf.grib','r') ! 打开数据源文件读取数据



 
  call codes_grib_new_from_file(datafile,igribdata,err)   ! 从数据源文件读取第一笔记录（记录号 igribdata) 
  call codes_get_size(igribdata,'values',size1)              ! 获取编号igribdata的记录中变量‘values’的大小

  allocate(v(size1))                                         ! 分配变量空间，用以存放读取的数据
  allocate(v1(size1))                                        ! 分配变量空间，用以存放用户创建的变量
  allocate(v2(size1))                                        ! 分配变量空间，用以存放用户创建的变量
 
  call codes_get(igribdata,'values',v)                       ! 从源文件读取数据，到变量v
 
  ! 构建新变量
  v=v*1000.0 ! different units for the output grib            
  v1=v
 

  do while (err/=CODES_END_OF_FILE)
 
    call codes_clone(igribsample,igribclone) ! clone sample before modifying it ! 复制示例文件的gribid
 
    ! 写入 关键字/变量 的值 （到加载的内存？）
    ! The given value is set for the key in the gribid message. In some cases the value can be an array rather than a scalar.
    ! As examples of array keys we have "values","pl", "pv" respectively the data values, 
    ! the list of number of points for each latitude in a reduced grid and the list of vertical levels.
    ! In these cases the value array must be allocated by the caller and their required dimension can be obtained with codes_get_size.
    !  为消息中的关键字写入给定值。在某些情况下，值可以是数组而不是标量。作为数组关键字的例子有“values”、“pl”、“pv”，分别表示是数据值、简化网格中每个纬度的点数列表和垂直层次列表。
    !  在这些情况下，需要使用函数codes_get_size来获得数组的大小。
    call codes_set(igribclone,'dataDate',date1)
    call codes_set(igribclone,'table2Version',table2Version)
    call codes_set(igribclone,'indicatorOfParameter',indicatorOfParameter) 
    call codes_set(igribclone,'stepType',stepType)
    call codes_set(igribclone,'startStep',startStep)
    call codes_set(igribclone,'endStep',endStep) 
    call codes_set(igribclone,'decimalPrecision',decimalPrecision) 
    call codes_set(igribclone,'values',v)


    ! 将更改写入到输出文件中 
    call codes_write(igribclone,outfile)
 
    ! 从源文件读取下一条记录
    call codes_grib_new_from_file(datafile,igribdata,err)
 
    if (err==0) then
      ! 读取变量'values'的值到v2中
      call codes_get(igribdata,'values',v2) 

      ! 输出单位乘以1000
      v2=v2*1000.0 ! different units for the output grib 
 
      ! 与初始时刻的变量值v1作差，计算起始步到终止步的累积量
      v=v2-v1 ! accumulation from startStep to endStep
 
      ! 时刻赋值，将本次步长的值赋值到v1中
      v1=v2 ! save previous step field
 
      ! 时间步长向后移动12h
      startStep=startStep+12
      endStep=endStep+12
 
    endif
 
  enddo
 
  ! 关闭示例的message
  call codes_release(igribsample)

  deallocate(v)
  deallocate(v1)
  deallocate(v2)
 
  ! 关闭输出文件
  call codes_close_file(outfile)
 
end program sample

 

参考自 training material

eccodes编码 使用eccodes 创建一个新的grib消息，包含集合平均的

创建grib消息有不同的选项

–  从其中一个输入grib场复制到要生成的新场。

– 使用默认示例目录中的样本（或模板）创建。参见“codes_info”。

– 使用私有示例目录中的样本创建。

• 第一个选项是最容易实现的。为了简单起见，我们建议您使用此选项。

•文件codes_create.f90（或grib_api_create.py)包含创建grib消息的代码框架。请试着添加创建一个grib消息所需的代码。按照代码中的说明操作。使用make编译代码，然后运行程序或运行Python脚本。

•注意！源代码中给出的更改GRIB头的参考文件已过时。请参阅下面http://apps.ecmwf.int/codes/grib/format/grib/local/ 关于本地定义 和 http://apps.ecmwf.int/codes/grib/format/mars/type/关于数据类型。

•现在将输入文件“eps”的链接更改为grib2文件（也可从~trx/ecCodes/data获得），然后再次运行该程序

eccodes_create_clone

! Copyright 2005-2019 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation 
! nor does it submit to any jurisdiction.
!
!
!  Description: How to create and use and index to access messages from a file
!               Creation of grib message, from a grib sample (or template).
!
!
!
!
program index
  use eccodes
  implicit none
  
  integer              :: iret
  integer,dimension(:),allocatable :: paramId,number
  character(len=20)    :: packingType
  character(len=5)     :: marsStream
  integer              :: paramIdSize,numberSize
  integer              :: i,j,is_missing
  integer              :: idx,igrib,igribout,count,numberOfValues
  integer              :: level,date,time,stepUnits,stepRange,startStep,endStep
  integer              :: localDefinitionNumber,Ni,Nj,PLPresent,N,nb_pl
  integer              :: generatingProcessIdentifier,numberOfForecastsInEnsemble
  integer              :: legBaseDate,legBaseTime,legNumber
  integer              :: bitmapPresent,bitsPerValue
  integer              :: outfile
  real (KIND=8),dimension(:), allocatable   :: values,result
  real (KIND=8)        :: latitudeOfFirstPoint,longitudeOfFirstPoint
  real (KIND=8)        :: latitudeOfLastPoint,longitudeOfLastPoint
  integer,dimension(:), allocatable   :: pl

  ! create an index for the file 'eps' file using the keys number and paramId
  ! use codes_index_create

  call codes_index_create(idx,'eps','paramId,number')

  ! get the number of distinct values of paramId in the index
  ! use codes_index_get_size

  call codes_index_get_size(idx,'paramId',paramIdSize)

  ! allocate the array to contain the list of distinct paramId

  allocate(paramId(paramIdSize))

  ! get the list of distinct paramId from the index
  ! use codes_index_get

  call codes_index_get(idx,'paramId',paramId)

  print*, 'grib contains ',paramIdSize, 'different parameters'

  ! get the number of distinct values of number in the index
  ! use codes_index_get_size

  call codes_index_get_size(idx,'number',numberSize)

  ! allocate the array to contain the list of distinct numbers

  allocate(number(numberSize))

  ! get the list of distinct numbers from the index
  ! use codes_index_get

  call codes_index_get(idx,'number',number)

  print*, 'grib contains ',numberSize, 'different EPS members'

  count=0

  ! select paramId 130 - Temperature
  ! use codes_index_select

  call codes_index_select(idx,'paramId','130')

  ! loop over the different ensemble members

  do j=1,numberSize ! loop on number

    ! select ensemble number=number(j)
    ! use codes_index_select

    call codes_index_select(idx,'number',number(j))

    ! get one grib message for the above values of the index
    ! use codes_new_from_index

    call codes_new_from_index(idx,igrib)
 
    ! for the first filed allocate array for values and result
    ! use codes_get_size

    if ( j == 1 ) then
      call codes_get_size(igrib,'values',numberOfValues)
      allocate(values(numberOfValues), stat = iret)
      if (iret /= 0)  STOP 'Failed to allocate values'
      allocate(result(numberOfValues), stat = iret)
      if (iret /= 0)  STOP 'Failed to allocate result'
    end if 

    ! get data values
    ! use codes_get

    call codes_get(igrib,"values", values)

    count = count + 1

    result=result+values

    ! release the grib message
    ! use grib_release
  
    if ( j /= numberSize ) then
       call codes_release(igrib)
    end if

  end do ! loop on number

  ! release the index
  ! use codes_index_release

  call codes_index_release(idx)

  print*,'We considered ',count,' members'

  result=result/count

  print*,'=============================================================================='
  print*, 'Stats for ensemble mean of T850'
  print*, 'Min: ', minval(result), 'Max: ', maxval(result), 'Avg: ', sum(result)/numberOfValues
  print*,'=============================================================================='

  ! take a sample grib message
  ! use codes_grib_new_from_samples:
  !   codes_grib_new_from_samples (gribid, samplename [, status])

  call codes_grib_new_from_samples(igribout, "reduced_gg_pl_grib1")

  ! open output grib file eps_mean.grib
  ! use codes_open_file:
  !   codes_open_file (ifile, filename, mode [, status])

  call codes_open_file(outfile, 'eps_mean.grib','w')

  ! change grib headers:
  !
  ! 1. product definition:
  !
  ! change grib headers
  ! See  http://apps.ecmwf.int/codes/grib/format/grib1/local/ (definition 1 or
  ! 30) and
  !      http://apps.ecmwf.int/codes/grib/format/mars/type/ for the datatype
  ! Keys to change are perturbationNumber and dataType.

  ! use codes_set:
  !   codes_set (gribid, key, value [, status])

  call codes_set(igribout,'paramId',130)
  call codes_get(igrib,'level',level)
  call codes_set(igribout,'level',level)
  call codes_get(igrib,'dataDate',date)
  call codes_set(igribout,'dataDate',date)
  call codes_get(igrib,'dataTime',time)
  call codes_set(igribout,'dataTime',time)
  call codes_get(igrib,'stepUnits',stepUnits)
  call codes_set(igribout,'stepUnits',stepUnits)
  call codes_get(igrib,'stepRange',stepRange)
  call codes_set(igribout,'stepRange',stepRange)
  call codes_get(igrib,'startStep',startStep)
  call codes_set(igribout,'startStep',startStep)
  call codes_get(igrib,'endStep',endStep)
  call codes_set(igribout,'endStep',endStep)
 
  call codes_set(igribout,'dataType',"em")
  call codes_get(igrib,'generatingProcessIdentifier',generatingProcessIdentifier)
  call codes_set(igribout,'generatingProcessIdentifier',generatingProcessIdentifier)
  call codes_get(igrib,'localDefinitionNumber',localDefinitionNumber)
  call codes_set(igribout,'localDefinitionNumber',localDefinitionNumber)
  call codes_get(igrib,'numberOfForecastsInEnsemble',numberOfForecastsInEnsemble)
  call codes_set(igribout,'numberOfForecastsInEnsemble',numberOfForecastsInEnsemble)
  call codes_get(igrib,'marsStream',marsStream)
  call codes_set(igribout,'marsStream',marsStream)

  call codes_get(igrib,'legBaseDate',legBaseDate)
  call codes_set(igribout,'legBaseDate',legBaseDate)
  call codes_get(igrib,'legBaseTime',legBaseTime)
  call codes_set(igribout,'legBaseTime',legBaseTime)
  call codes_get(igrib,'legNumber',legNumber)
  call codes_set(igribout,'legNumber',legNumber)

  ! 2. grid definition:
  !
  ! see https://software.ecmwf.int/wiki/display/GRIB/Gaussian+grids

  is_missing=0;
  call codes_is_missing(igrib,'Ni',is_missing);
  if ( is_missing /= 1 ) then
     call codes_get(igrib,'Ni',Ni) 
     call codes_set(igribout,'Ni',Ni)
  endif     

  is_missing=0;
  call codes_is_missing(igrib,'Nj',is_missing);
  if ( is_missing /= 1 ) then
     call codes_get(igrib,'Nj',Nj) 
     call codes_set(igribout,'Nj',Nj)
  endif     

  call codes_get(igrib,'latitudeOfFirstGridPointInDegrees',latitudeOfFirstPoint)
  call codes_set(igribout,'latitudeOfFirstGridPointInDegrees',latitudeOfFirstPoint)
  call codes_get(igrib,'longitudeOfFirstGridPointInDegrees',longitudeOfFirstPoint)
  call codes_set(igribout,'longitudeOfFirstGridPointInDegrees',longitudeOfFirstPoint)
  call codes_get(igrib,'latitudeOfLastGridPointInDegrees',latitudeOfLastPoint)
  call codes_set(igribout,'latitudeOfLastGridPointInDegrees',latitudeOfLastPoint)
  call codes_get(igrib,'longitudeOfLastGridPointInDegrees',longitudeOfLastPoint)
  call codes_set(igribout,'longitudeOfLastGridPointInDegrees',longitudeOfLastPoint)
  call codes_get(igrib,'N',N)
  call codes_set(igribout,'N',N)
  call codes_get(igrib,'PLPresent',PLPresent)
  if (PLPresent == 1) then
     call codes_get_size(igrib,'pl',nb_pl)
     allocate(pl(nb_pl))
     call codes_get(igrib,'pl',pl)
     call codes_set(igribout,'pl',pl)
     deallocate(pl)
  else
     print*, "There is no PL values in your GRIB message!"
  end if
  
  !  3. bitmap definition:

  call codes_get(igrib,'bitmapPresent',bitmapPresent)
  call codes_set(igribout,'bitmapPresent',bitmapPresent)

  ! 4. data:

  call codes_get(igrib,'bitsPerValue',bitsPerValue)
  call codes_set(igribout,'bitsPerValue',bitsPerValue)
  call codes_get(igrib,'packingType',packingType)
  call codes_set(igribout,'packingType',packingType)

  call codes_set(igribout,'values',result)

  ! write message out
  ! use codes_write:
  !   codes_write (gribid, ifile [, status])

  call codes_write(igribout,outfile)

  ! release the output grib message
  ! use grib_release

  call codes_release(igribout)


end program index

 

eccoded_create_sample

! Copyright 2005-2018 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation 
! nor does it submit to any jurisdiction.
!
!
!  Description: How to create and use and index to access messages from a file
!               Creation of grib message, from a grib sample (or template).
!
!
!
!
program index
  use eccodes
  implicit none
  
  integer              :: iret
  integer,dimension(:),allocatable :: paramId,number
  character(len=20)    :: packingType
  character(len=5)     :: marsStream
  integer              :: paramIdSize,numberSize
  integer              :: i,j,is_missing
  integer              :: idx,igrib,igribout,count,numberOfValues
  integer              :: level,date,time,stepUnits,stepRange,startStep,endStep
  integer              :: localDefinitionNumber,Ni,Nj,PLPresent,N,nb_pl
  integer              :: generatingProcessIdentifier,numberOfForecastsInEnsemble
  integer              :: legBaseDate,legBaseTime,legNumber
  integer              :: bitmapPresent,bitsPerValue
  integer              :: outfile
  real (KIND=8),dimension(:), allocatable   :: values,result
  real (KIND=8)        :: latitudeOfFirstPoint,longitudeOfFirstPoint
  real (KIND=8)        :: latitudeOfLastPoint,longitudeOfLastPoint
  integer,dimension(:), allocatable   :: pl

  ! create an index for the file 'eps' file using the keys number and paramId
  ! use codes_index_create

  call codes_index_create(idx,'eps','paramId,number')

  ! get the number of distinct values of paramId in the index
  ! use codes_index_get_size

  call codes_index_get_size(idx,'paramId',paramIdSize)

  ! allocate the array to contain the list of distinct paramId

  allocate(paramId(paramIdSize))

  ! get the list of distinct paramId from the index
  ! use codes_index_get

  call codes_index_get(idx,'paramId',paramId)

  print*, 'grib contains ',paramIdSize, 'different parameters'

  ! get the number of distinct values of number in the index
  ! use codes_index_get_size

  call codes_index_get_size(idx,'number',numberSize)

  ! allocate the array to contain the list of distinct numbers

  allocate(number(numberSize))

  ! get the list of distinct numbers from the index
  ! use codes_index_get

  call codes_index_get(idx,'number',number)

  print*, 'grib contains ',numberSize, 'different EPS members'

  count=0

  ! select paramId 130 - Temperature
  ! use codes_index_select

  call codes_index_select(idx,'paramId','130')

  ! loop over the different ensemble members

  do j=1,numberSize ! loop on number

    ! select ensemble number=number(j)
    ! use codes_index_select

    call codes_index_select(idx,'number',number(j))

    ! get one grib message for the above values of the index
    ! use codes_new_from_index

    call codes_new_from_index(idx,igrib)
 
    ! for the first filed allocate array for values and result
    ! use codes_get_size

    if ( j == 1 ) then
      call codes_get_size(igrib,'values',numberOfValues)
      allocate(values(numberOfValues), stat = iret)
      if (iret /= 0)  STOP 'Failed to allocate values'
      allocate(result(numberOfValues), stat = iret)
      if (iret /= 0)  STOP 'Failed to allocate result'
    end if 

    ! get data values
    ! use codes_get

    call codes_get(igrib,"values", values)

    count = count + 1

    result=result+values

    ! release the grib message
    ! use grib_release
  
    if ( j /= numberSize ) then
       call codes_release(igrib)
    end if

  end do ! loop on number

  ! release the index
  ! use codes_index_release

  call codes_index_release(idx)

  print*,'We considered ',count,' members'

  result=result/count

  print*,'=============================================================================='
  print*, 'Stats for ensemble mean of T850'
  print*, 'Min: ', minval(result), 'Max: ', maxval(result), 'Avg: ', sum(result)/numberOfValues
  print*,'=============================================================================='

  ! take a sample grib message
  ! use codes_grib_new_from_samples:
  !   codes_grib_new_from_samples (gribid, samplename [, status])

  call codes_grib_new_from_samples(igribout, "reduced_gg_pl_grib1")

  ! open output grib file eps_mean.grib
  ! use codes_open_file:
  !   codes_open_file (ifile, filename, mode [, status])

  call codes_open_file(outfile, 'eps_mean.grib','w')

  ! change grib headers:
  !
  ! 1. product definition:
  !
  ! change grib headers
  ! See  http://apps.ecmwf.int/codes/grib/format/grib1/local/ (definition 1 or
  ! 30) and
  !      http://apps.ecmwf.int/codes/grib/format/mars/type/ for the datatype
  ! Keys to change are perturbationNumber and dataType.

  ! use codes_set:
  !   codes_set (gribid, key, value [, status])

  call codes_set(igribout,'paramId',130)
  call codes_get(igrib,'level',level)
  call codes_set(igribout,'level',level)
  call codes_get(igrib,'dataDate',date)
  call codes_set(igribout,'dataDate',date)
  call codes_get(igrib,'dataTime',time)
  call codes_set(igribout,'dataTime',time)
  call codes_get(igrib,'stepUnits',stepUnits)
  call codes_set(igribout,'stepUnits',stepUnits)
  call codes_get(igrib,'stepRange',stepRange)
  call codes_set(igribout,'stepRange',stepRange)
  call codes_get(igrib,'startStep',startStep)
  call codes_set(igribout,'startStep',startStep)
  call codes_get(igrib,'endStep',endStep)
  call codes_set(igribout,'endStep',endStep)
 
  call codes_set(igribout,'dataType',"em")
  call codes_get(igrib,'generatingProcessIdentifier',generatingProcessIdentifier)
  call codes_set(igribout,'generatingProcessIdentifier',generatingProcessIdentifier)
  call codes_get(igrib,'localDefinitionNumber',localDefinitionNumber)
  call codes_set(igribout,'localDefinitionNumber',localDefinitionNumber)
  call codes_get(igrib,'numberOfForecastsInEnsemble',numberOfForecastsInEnsemble)
  call codes_set(igribout,'numberOfForecastsInEnsemble',numberOfForecastsInEnsemble)
  call codes_get(igrib,'marsStream',marsStream)
  call codes_set(igribout,'marsStream',marsStream)

  call codes_get(igrib,'legBaseDate',legBaseDate)
  call codes_set(igribout,'legBaseDate',legBaseDate)
  call codes_get(igrib,'legBaseTime',legBaseTime)
  call codes_set(igribout,'legBaseTime',legBaseTime)
  call codes_get(igrib,'legNumber',legNumber)
  call codes_set(igribout,'legNumber',legNumber)

  ! 2. grid definition:
  !
  ! see https://software.ecmwf.int/wiki/display/GRIB/Gaussian+grids

  is_missing=0;
  call codes_is_missing(igrib,'Ni',is_missing);
  if ( is_missing /= 1 ) then
     call codes_get(igrib,'Ni',Ni) 
     call codes_set(igribout,'Ni',Ni)
  endif     

  is_missing=0;
  call codes_is_missing(igrib,'Nj',is_missing);
  if ( is_missing /= 1 ) then
     call codes_get(igrib,'Nj',Nj) 
     call codes_set(igribout,'Nj',Nj)
  endif     

  call codes_get(igrib,'latitudeOfFirstGridPointInDegrees',latitudeOfFirstPoint)
  call codes_set(igribout,'latitudeOfFirstGridPointInDegrees',latitudeOfFirstPoint)
  call codes_get(igrib,'longitudeOfFirstGridPointInDegrees',longitudeOfFirstPoint)
  call codes_set(igribout,'longitudeOfFirstGridPointInDegrees',longitudeOfFirstPoint)
  call codes_get(igrib,'latitudeOfLastGridPointInDegrees',latitudeOfLastPoint)
  call codes_set(igribout,'latitudeOfLastGridPointInDegrees',latitudeOfLastPoint)
  call codes_get(igrib,'longitudeOfLastGridPointInDegrees',longitudeOfLastPoint)
  call codes_set(igribout,'longitudeOfLastGridPointInDegrees',longitudeOfLastPoint)
  call codes_get(igrib,'N',N)
  call codes_set(igribout,'N',N)
  call codes_get(igrib,'PLPresent',PLPresent)
  if (PLPresent == 1) then
     call codes_get_size(igrib,'pl',nb_pl)
     allocate(pl(nb_pl))
     call codes_get(igrib,'pl',pl)
     call codes_set(igribout,'pl',pl)
     deallocate(pl)
  else
     print*, "There is no PL values in your GRIB message!"
  end if
  
  !  3. bitmap definition:

  call codes_get(igrib,'bitmapPresent',bitmapPresent)
  call codes_set(igribout,'bitmapPresent',bitmapPresent)

  ! 4. data:

  call codes_get(igrib,'bitsPerValue',bitsPerValue)
  call codes_set(igribout,'bitsPerValue',bitsPerValue)
  call codes_get(igrib,'packingType',packingType)
  call codes_set(igribout,'packingType',packingType)

  call codes_set(igribout,'values',result)

  ! write message out
  ! use codes_write:
  !   codes_write (gribid, ifile [, status])

  call codes_write(igribout,outfile)

  ! release the output grib message
  ! use grib_release

  call codes_release(igribout)


end program index