00001 ! Copyright 2005-2007 ECMWF 00002 ! 00003 ! Licensed under the GNU Lesser General Public License which 00004 ! incorporates the terms and conditions of version 3 of the GNU 00005 ! General Public License. 00006 ! See LICENSE and gpl-3.0.txt for details. 00007 ! 00008 ! 00009 ! Description: how to create a new GRIB message from a template. 00010 ! 00011 ! 00012 ! 00013 ! Author: Anne Fouilloux 00014 ! 00015 ! 00016 ! 00017 program template 00018 use grib_api 00019 implicit none 00020 integer, parameter :: NbPressureLevel = 10 00021 integer :: err 00022 integer :: nx, ny, i 00023 integer :: outfile 00024 integer :: igrib 00025 integer :: date, level 00026 real :: missingValue 00027 integer :: indicatorOfParameter 00028 double precision, dimension(:,:), allocatable :: field2D 00029 integer, dimension(NbPressureLevel) :: levels=(/1,2,10,50,70,100,200,500,850,1000/) 00030 00031 ! Initialisation 00032 00033 date = 20080104 00034 indicatorOfParameter = 130 ! temperature in table 128 used in the template 00035 00036 ! a new grib message is loaded from an existing template 00037 ! templates are searched in a default template path (use grib_info 00038 ! to see where it is. The default template path can be changed by 00039 ! setting the environment variable GRIB_TEMPLATES_PATH 00040 call grib_new_from_template(igrib, "GRIB1") 00041 00042 ! Open a new GRIB file to write new GRIB messages 00043 call grib_open_file(outfile, 'out.grib1','w') 00044 00045 call grib_get(igrib,"numberOfPointsAlongAParallel", nx) 00046 00047 call grib_get(igrib,"numberOfPointsAlongAMeridian",ny) 00048 00049 call grib_get(igrib,'lev',level) 00050 00051 allocate(field2D(nx,ny),stat=err) 00052 00053 if (err .ne. 0) then 00054 print*, 'Failed to allocate ', nx*ny, ' values' 00055 STOP 00056 end if 00057 call generate_field(field2D, missingValue) 00058 00059 ! Same date and same parameters for all the levels 00060 call grib_set(igrib,"dataDate", date) 00061 call grib_set(igrib,"indicatorOfParameter", indicatorOfParameter) 00062 00063 ! Write new GRIBs from the same message for different pressure levels 00064 do i=1, size(levels) 00065 call grib_set(igrib,'lev',levels(i)) 00066 00067 ! compute the field for each pressure level 00068 call generate_field(field2D, missingValue) 00069 00070 ! set the value of missingValue 00071 call grib_set(igrib, 'missingValue',missingValue) 00072 00073 call grib_set(igrib,"bitmapPresent",1) 00074 00075 ! use pack to create 1D values 00076 call grib_set(igrib,'values',pack(field2D, mask=.true.)) 00077 00078 ! write the new message to a file 00079 call grib_write(igrib,outfile) 00080 end do 00081 00082 call grib_release(igrib) 00083 00084 call grib_close_file(outfile) 00085 deallocate(field2D) 00086 00087 contains 00088 !====================================== 00089 subroutine generate_field(gfield2D, missingValue) 00090 double precision, dimension(:,:) :: gfield2D 00091 real, intent(out) :: missingValue 00092 00093 integer, parameter :: borders = 5 00094 integer :: i 00095 integer, dimension(2) :: shape_gfield2D 00096 00097 shape_gfield2D = shape(gfield2D) 00098 00099 call random_number(gfield2D) 00100 ! define the value of missingValue 00101 missingValue = 99999 00102 ! just keep the boundaries for gfield2D 00103 gfield2D(borders+1:shape_gfield2D(1)-borders,borders+1:shape_gfield2D(2)-borders) = missingValue 00104 00105 end subroutine generate_field 00106 !====================================== 00107 00108 end program template