precision_fortran.F

precision_fortran.F How to control precision when coding a grib field.

00001 C Copyright 2005-2016 ECMWF
00002 C This software is licensed under the terms of the Apache Licence Version 2.0
00003 C which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
00004 C 
00005 C In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
00006 C virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
00007 C
00008 C
00009 C  Fortran 77 Implementation: precision
00010 C
00011 C  Description: how to control decimal precision when packing fields.
00012 C
00013 C
00014 C  Author: Enrico Fucile
00015 C
00016 C
00017 C
00018       program precision
00019       implicit none
00020       integer maxNumberOfValues
00021       parameter (maxNumberOfValues=10000)
00022       include 'grib_api_f77.h'
00023       integer*4 size
00024       integer infile,outfile
00025       integer igrib
00026       real*8 values1(maxNumberOfValues)
00027       real*8 values2(maxNumberOfValues)
00028       real*8 maxa,a,maxv,minv,maxr,r
00029       integer*4 decimalPrecision,bitsPerValue1,bitsPerValue2
00030       integer i
00031 
00032       call grib_check(grib_open_file(infile
00033      X,'../../data/regular_latlon_surface.grib1','r'))
00034 
00035       call grib_check(grib_open_file(outfile
00036      X,'../../data/regular_latlon_surface_prec.grib1','w'))
00037 
00038 C     a new grib message is loaded from file
00039 C     igrib is the grib id to be used in subsequent calls
00040       call grib_check(grib_new_from_file(infile,igrib))
00041 
00042 C     bitsPerValue before changing the packing parameters
00043       call grib_check(grib_get_int(igrib,'bitsPerValue',bitsPerValue1))
00044 
00045 C     get the size of the values array
00046       call grib_check(grib_get_size(igrib,"values",size))
00047 
00048 C     get data values before changing the packing parameters*/
00049       call grib_check(grib_get_real8_array(igrib,"values",values1,size))
00050 
00051 C     setting decimal precision=2 means that 2 decimal digits
00052 C     are preserved when packing.
00053       decimalPrecision=2
00054       call grib_check(grib_set_int(igrib,"changeDecimalPrecision"
00055      X,decimalPrecision))
00056 
00057 C     bitsPerValue after changing the packing parameters
00058       call grib_check(grib_get_int(igrib,"bitsPerValue",bitsPerValue2))
00059 
00060 C     get data values after changing the packing parameters
00061       call grib_check(grib_get_real8_array(igrib,"values",values2,size))
00062 
00063 C     computing error
00064       maxa=0
00065       maxr=0
00066       maxv=values2(1)
00067       minv=maxv
00068       do i=1,size
00069         a=abs(values2(i)-values1(i))
00070         if ( values2(i) .gt. maxv ) maxv=values2(i)
00071         if ( values2(i) .lt. maxv ) minv=values2(i)
00072         if ( values2(i) .ne. 0 ) then
00073          r=abs((values2(i)-values1(i))/values2(i))
00074         endif
00075         if ( a .gt. maxa ) maxa=a
00076         if ( r .gt. maxr ) maxr=r
00077       enddo
00078       write(*,*) "max absolute error = ",maxa
00079       write(*,*) "max relative error = ",maxr
00080       write(*,*) "min value = ",minv
00081       write(*,*) "max value = ",maxv
00082 
00083       write(*,*) "old number of bits per value=",bitsPerValue1
00084       write(*,*) "new number of bits per value=",bitsPerValue2
00085 
00086 C     write modified message to a file
00087       call grib_check(grib_write(igrib,outfile))
00088 
00089       call grib_check(grib_release(igrib))
00090 
00091       call grib_check(grib_close_file(infile))
00092 
00093       call grib_check(grib_close_file(outfile))
00094 
00095       end
00096 

Generated on Tue Sep 22 15:18:21 2009 for grib_api by  doxygen 1.5.3