You are viewing an old version of this page. View the current version.

Compare with Current View Page History

« Previous Version 3 Next »

Description

This example shows: how to set pv values in a GRIB message

Source code

grib_set_pv.f90
! Copyright 2005-2016 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 set pv values in a GRIB message
!
!
program grib_set_pv
  use eccodes
  implicit none
  integer                         :: numberOfLevels
  integer                         :: numberOfCoefficients
  integer                         :: outfile, igrib
  integer                         :: i, ios
  real, dimension(:),allocatable  :: pv
  
  numberOfLevels=60
  numberOfCoefficients=2*(numberOfLevels+1)

  allocate(pv(numberOfCoefficients))

  ! read the model level coefficients from file
  open( unit=1, file="../../data/60_model_levels", &
                form="formatted",action="read")

  do i=1,numberOfCoefficients,2
     read(unit=1,fmt=*, iostat=ios) pv(i), pv(i+1)
     if (ios /= 0) then
        print *, "I/O error: ",ios
        exit
     end if
  end do
  
  ! print coefficients
  !do i=1,numberOfCoefficients,2
  !  print *,"  a=",pv(i)," b=",pv(i+1)
  !end do

  close(unit=1)

  call codes_open_file(outfile, 'out.pv.grib1','w')
  
  !     a new grib message is loaded from file
  !     igrib is the grib id to be used in subsequent calls
  call codes_new_from_samples(igrib, "reduced_gg_sfc_grib1")

  !     set levtype to ml (model level)
  call codes_set(igrib,'typeOfLevel','hybrid')

  !     set level 
  call codes_set(igrib,'level',2)

  !     set PVPresent as an integer 
  call codes_set(igrib,'PVPresent',1)
  
  call codes_set(igrib,'pv',pv)
  
  !     write modified message to a file
  call codes_write(igrib,outfile)
  
  !  FREE MEMORY
  call codes_release(igrib)
  deallocate(pv)

  call codes_close_file(outfile)
  
end program grib_set_pv
grib_set_pv.py
#
# Copyright 2005-2016 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.
#

import traceback
import sys

from eccodes import *

VERBOSE = 1  # verbose error reporting


def example():
    """
    Encoding of the pv coefficients.
    """
    # read the coefficients from file
    pv = []
    for line in open('../../data/60_model_levels'):
        pv.extend([float(x) for x in line.strip().split('\t')])

    numberOfLevels = 60
    numberOfCoefficients = 2 * (numberOfLevels + 1)
    assert(len(pv) == numberOfCoefficients)

    fout = open('out.pv.grib1', 'w')
    gid = codes_new_from_samples('reduced_gg_sfc_grib1')

    codes_set(gid, 'typeOfLevel', 'hybrid')
    codes_set(gid, 'level', 2)
    codes_set(gid, 'PVPresent', 1)
    codes_set_array(gid, 'pv', pv)

    codes_write(gid, fout)

    codes_release(gid)
    fout.close()


def main():
    try:
        example()
    except CodesInternalError, err:
        if VERBOSE:
            traceback.print_exc(file=sys.stderr)
        else:
            print >>sys.stderr, err.msg

        return 1

if __name__ == '__main__':
    sys.exit(main())
grib_set_pv.c
/*
 * Copyright 2005-2016 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.
 */

/*
 * C Implementation: grib_set_pv
 *
 * Description: how to set pv (vertical coordinate parameters)
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "eccodes.h"

void usage(const char* prog) {
    fprintf(stderr, "usage: %s in out\n",prog);
    exit(1);
}

int main(int argc, char** argv)
{
    int err = 0;
    long NV = 0;
    size_t size=0;
    double pv[4]={1,2,3,4};
    size_t pvsize=4;

    FILE* in = NULL;
    char* infile = NULL;
    FILE* out = NULL;
    char* outfile = NULL;
    codes_handle *h = NULL;
    const void* buffer = NULL;

    if (argc != 3) usage(argv[0]);
    infile = argv[1];
    outfile= argv[2];

    in = fopen(infile, "r");
    if(!in) {
        fprintf(stderr, "ERROR: unable to open input file %s\n",infile);
        return 1;
    }

    out = fopen(outfile, "w");
    if(!out) {
        fprintf(stderr, "ERROR: unable to open output file %s\n",outfile);
        fclose(in);
        return 1;
    }

    /* create a new handle from a message in a file */
    h = codes_handle_new_from_file(0, in, PRODUCT_GRIB, &err);
    if (h == NULL) {
        fprintf(stderr, "Error: unable to create handle from file %s\n",infile);
    }

    CODES_CHECK(codes_set_long(h,"PVPresent", 1),0);

    CODES_CHECK(codes_set_double_array(h, "pv", pv, pvsize),0);

    /* Once we set the pv array, the NV key should be also set */
    CODES_CHECK(codes_get_long(h, "NV", &NV),0);
    assert( NV == pvsize );

    /* get the coded message in a buffer */
    CODES_CHECK(codes_get_message(h, &buffer, &size),0);

    /* write the buffer in a file*/
    if(fwrite(buffer, 1, size, out) != size)
    {
        perror(argv[1]);
        exit(1);
    }

    codes_handle_delete(h);
    fclose(in);
    fclose(out);

    return 0;
}
  • No labels