nearest.f90

How to find the nearest grid points.

00001 ! Copyright 2005-2016 ECMWF
00002 ! This software is licensed under the terms of the Apache Licence Version 2.0
00003 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
00004 ! 
00005 ! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
00006 ! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
00007 !
00008 !
00009 !  Description: how to use grib_find_nearest and grib_get_element 
00010 !
00011 !
00012 !  Author: Enrico Fucile 
00013 !
00014 !
00015 !
00016 program find
00017   use grib_api
00018   implicit none
00019   integer                                      :: npoints
00020   integer                                      :: infile
00021   integer                                      :: igrib, ios, i
00022   real(8), dimension(:), allocatable  :: lats, lons
00023   real(8), dimension(:), allocatable  :: nearest_lats, nearest_lons
00024   real(8), dimension(:), allocatable  :: distances, values, lsm_values
00025   integer(kind=kindOfInt), dimension(:), allocatable  :: indexes
00026   real(kind=8)                        :: value
00027 
00028 ! initialization
00029   open( unit=1, file="../../data/list_points",form="formatted",action="read")
00030   read(unit=1,fmt=*) npoints
00031   allocate(lats(npoints))
00032   allocate(lons(npoints))
00033   allocate(nearest_lats(npoints))
00034   allocate(nearest_lons(npoints))
00035   allocate(distances(npoints))
00036   allocate(lsm_values(npoints))
00037   allocate(values(npoints))
00038   allocate(indexes(npoints))
00039   do i=1,npoints
00040      read(unit=1,fmt=*, iostat=ios) lats(i), lons(i)
00041      if (ios /= 0) then
00042         npoints = i - 1
00043         exit
00044      end if
00045   end do
00046   close(unit=1)
00047   call grib_open_file(infile, &
00048        '../../data/reduced_gaussian_lsm.grib1','r')
00049   
00050   !     a new grib message is loaded from file
00051   !     igrib is the grib id to be used in subsequent calls
00052   call grib_new_from_file(infile,igrib)
00053   
00054 
00055   call grib_find_nearest(igrib, .true., lats, lons, nearest_lats, nearest_lons,lsm_values, distances, indexes)
00056   call grib_release(igrib)
00057   
00058   call grib_close_file(infile)
00059 
00060 ! will apply it to another GRIB
00061   call grib_open_file(infile, &
00062        '../../data/reduced_gaussian_pressure_level.grib1','r')
00063   call grib_new_from_file(infile,igrib)
00064 
00065   call grib_get_element(igrib,"values", indexes, values)
00066   call grib_release(igrib)
00067   call grib_close_file(infile)
00068 
00069   do i=1, npoints
00070      print*,lats(i), lons(i), nearest_lats(i), nearest_lons(i), distances(i), lsm_values(i), values(i)
00071   end do
00072 
00073   deallocate(lats)
00074   deallocate(lons)
00075   deallocate(nearest_lats)
00076   deallocate(nearest_lons)
00077   deallocate(distances)
00078   deallocate(lsm_values)
00079   deallocate(values)
00080   deallocate(indexes)
00081 
00082 end program find

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