1! (C) Copyright 2005- ECMWF.
2!
3! This software is licensed under the terms of the Apache Licence Version 2.0
4! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5!
6! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
7! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
8!
9!
10!  Description: How to set pv values in a GRIB message
11!
12!
13program grib_set_pv
14   use eccodes
15   implicit none
16   integer                         :: numberOfLevels
17   integer                         :: numberOfCoefficients
18   integer                         :: outfile, igrib
19   integer                         :: i, ios
20   real, dimension(:), allocatable  :: pv
21
22   numberOfLevels = 60
23   numberOfCoefficients = 2*(numberOfLevels + 1)
24
25   allocate (pv(numberOfCoefficients))
26
27   ! read the model level coefficients from file
28   open (unit=1, file="../../data/60_model_levels", &
29         form="formatted", action="read")
30
31   do i = 1, numberOfCoefficients, 2
32      read (unit=1, fmt=*, iostat=ios) pv(i), pv(i + 1)
33      if (ios /= 0) then
34         print *, "I/O error: ", ios
35         exit
36      end if
37   end do
38
39   ! print coefficients
40   !do i=1,numberOfCoefficients,2
41   !  print *,"  a=",pv(i)," b=",pv(i+1)
42   !end do
43
44   close (unit=1)
45
46   call codes_open_file(outfile, 'out.pv.grib1', 'w')
47
48   ! A new grib message is loaded from file
49   ! igrib is the grib id to be used in subsequent calls
50   call codes_grib_new_from_samples(igrib, "reduced_gg_sfc_grib1")
51
52   ! set levtype to ml (model level)
53   call codes_set(igrib, 'typeOfLevel', 'hybrid')
54
55   ! set level
56   call codes_set(igrib, 'level', 2)
57
58   ! set PVPresent as an integer
59   call codes_set(igrib, 'PVPresent', 1)
60
61   call codes_set(igrib, 'pv', pv)
62
63   ! write modified message to a file
64   call codes_write(igrib, outfile)
65
66   ! Free memory
67   call codes_release(igrib)
68   deallocate (pv)
69
70   call codes_close_file(outfile)
71
72end program grib_set_pv
73