This Fortran90 example uses GRIB_API to read a GRIB edition 1 data file containing fields in Reduced Gaussian grid and interpolate these to a 0.5/0.5 degree LatLon grid.
Fortran 90 code
! ! Copyright 2015 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. ! ! Unless required by applicable law or agreed to in writing, software ! distributed under the License is distributed on an "AS IS" BASIS, ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ! ! 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. ! PROGRAM INTF_EXAMPLE_WITH_ARRAY ! USE grib_api IMPLICIT NONE INTEGER, DIMENSION(4) :: INTV REAL, DIMENSION(4) :: REALV CHARACTER(Len=20), DIMENSION(4) :: CHARV ! INTEGER (KIND=4), DIMENSION(1) :: INGRIB(1), NEWFLD(1) ! INTEGER :: IUNIT1, IUNIT2 INTEGER :: IGRIB INTEGER :: IRET, NARGS INTEGER (KIND=KINDOFINT) :: INLEN, NEWLEN INTEGER (KIND=KINDOFSIZE) :: INSIZE, NEWSIZE INTEGER :: J INTEGER :: NRESOL, NPL INTEGER :: Ni, Nj INTEGER :: IBPV INTEGER :: len_in, len_out INTEGER :: IPARAM, ITABLE INTEGER, DIMENSION(:), ALLOCATABLE :: NDGL REAL, DIMENSION(:), ALLOCATABLE :: VALUES_IN, VALUES_OUT REAL :: DLAT, DLON REAL :: VALMISS ! ! Externals INTEGER :: INTIN, INTOUT, INTF, IARGC CHARACTER(Len=128) :: INFILE, OUTFILE CHARACTER(Len=128), DIMENSION(20) :: CARG CHARACTER(Len=128) :: MESSAGE = ' ' CHARACTER(Len=32) :: CGTYPE = ' ', CINGRIDT = ' ' ! ! ********************************************************************** ! Default values IBPV = 16 ! bitsPerValue / accuracy NRESOL = 0 ! Gaussian resolution DLAT = 0.5 ! Latitude increment (in degrees) DLON = 0.5 ! Longitude increment (in degrees) VALMISS = -99999.9 ! Missing value indicator ! ! Pick up file names from command line. ! NARGS = 0 NARGS = IARGC() IF ( NARGS .NE. 4) THEN PRINT *, 'NARGS is incorrect' CALL USAGE STOP END IF DO J = 1, NARGS CALL GETARG(J,CARG(J)) END DO DO J = 1, NARGS,2 IF (CARG(J) == '-i') THEN INFILE=TRIM(CARG(J+1)) ELSEIF (CARG(J) == '-o') THEN OUTFILE=TRIM(CARG(J+1)) ELSE CALL USAGE STOP END IF END DO ! WRITE(*,*) '*** Open input and output files' ! CALL GRIB_OPEN_FILE (IUNIT1, INFILE, 'r', IRET) IF ( IRET /= GRIB_SUCCESS) STOP ' GRIB_OPEN of INFILE failed' ! CALL GRIB_OPEN_FILE (IUNIT2, OUTFILE, 'w', IRET) IF ( IRET /= GRIB_SUCCESS) STOP ' GRIB_OPEN of OUTFILE failed' ! WRITE(*,*) '*** Start of loop through input GRIB-coded fields' ! CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) LOOP: DO WHILE (IRET /= GRIB_END_OF_FILE) IF (IRET /= GRIB_SUCCESS) THEN PRINT *, ' GRIB_NEW_FROM_FILE failed with error code = ', IRET CALL GRIB_GET_ERROR_STRING(IRET,MESSAGE) PRINT *, TRIM(MESSAGE) PRINT *, 'Skipping to next message ...' CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) CYCLE ENDIF ! WRITE(*,*) '*** Unpack the values in the input GRIB' ! CALL GRIB_GET_SIZE(IGRIB,'values',len_in) WRITE(*,*) ' Allocating array for input field of size=',len_in ALLOCATE(VALUES_IN(len_in)) CALL GRIB_SET(IGRIB,'missingValue',VALMISS) CALL GRIB_GET(IGRIB,'values',VALUES_IN) ! ! Define the input field via calls to INTIN ! CHARV(1) = 'unpacked' IRET = INTIN('form',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTIN failed to set format' STOP ENDIF ! Get parameter information (assumes GRIB edition 1 here !) CALL GRIB_GET(IGRIB,'indicatorOfParameter',IPARAM) INTV(1) = IPARAM IRET = INTIN('parameter',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTIN failed to set parameter' STOP ENDIF CALL GRIB_GET(IGRIB,'table2Version',ITABLE) INTV(1) = ITABLE IRET = INTIN('table',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTIN failed to set table' STOP ENDIF ! Check it's a reduced Gaussian grid and set resolution CALL GRIB_GET(IGRIB,'gridType',CINGRIDT) IF (TRIM(CINGRIDT) == 'reduced_gg') THEN CALL GRIB_GET(IGRIB,'N',NRESOL) INTV(1) = NRESOL WRITE(*,*) ' Input resolution = ', NRESOL IRET = INTIN('reduced',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTIN failed to set reduced Gaussian resolution' STOP ENDIF ELSE WRITE(*,*) 'Input grid type ', TRIM(CINGRIDT), ' not reduced_gg' CYCLE ENDIF ! ! Set the missing value ! CHARV(1) = 'yes' REALV(1) = VALMISS IRET = INTIN('missingvalue',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTIN failed to set missingValue' STOP ENDIF ! ! Define the output field with calls to INTOUT ! CHARV(1) = 'unpacked' IRET = INTOUT('form',INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTOUT failed to set format' STOP ENDIF ! WRITE(*,*) '*** Define the grid interval for the new field(s)' REALV(1) = DLAT REALV(2) = DLON IRET = INTOUT('grid', INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTOUT failed to set grid resolution' STOP ENDIF WRITE(*,*) ' GRID TYPE = lat-lon' WRITE(*,*) ' RESOL = ', REALV(1), 'x', REALV(2) ! Calculate the number of points (assumes global grid) Ni = INT(360.0/REALV(2)) Nj = INT(180.0/REALV(1)) + 1 len_out = Ni*Nj WRITE(*,*) ' Allocating output array of size = ', len_out ALLOCATE(VALUES_OUT(LEN_OUT)) ! ! Define the packing accuracy for the new field(s). ! INTV(1) = IBPV IRET = INTOUT('accuracy', INTV, REALV, CHARV) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTOUT failed to set accuracy' STOP ENDIF WRITE(*,*) ' ACCURACY = ', IBPV ! WRITE(*,*) '*** Interpolate' INLEN =1 NEWLEN = 1 IRET = INTF(INGRIB,len_in,VALUES_IN,NEWFLD,len_out,VALUES_OUT) IF ( IRET /= 0 ) THEN WRITE(*,*) ' INTF failed for product ' WRITE(*,*) ' Skipping to next message ...' CALL GRIB_RELEASE(IGRIB) DEALLOCATE(VALUES_IN,VALUES_OUT) CYCLE ENDIF ! WRITE(*,*) '*** Set-up and write output GRIB' ! CALL GRIB_SET(IGRIB,'gridType','regular_ll') CALL GRIB_SET(IGRIB,'latitudeOfFirstGridPointInDegrees',90.0) CALL GRIB_SET(IGRIB,'latitudeOfLastGridPointInDegrees',-90.0) CALL GRIB_SET(IGRIB,'longitudeOfFirstGridPointInDegrees',0.0) CALL GRIB_SET(IGRIB,'longitudeOfLastGridPointInDegrees',360.0 - DLON) CALL GRIB_SET(IGRIB,'iDirectionIncrementInDegrees',DLON) CALL GRIB_SET(IGRIB,'jDirectionIncrementInDegrees',DLON) CALL GRIB_SET(IGRIB,'bitsPerValue',IBPV) CALL GRIB_SET(IGRIB,'Ni',Ni) CALL GRIB_SET(IGRIB,'Nj',Nj) CALL GRIB_SET(IGRIB,'values',VALUES_OUT) CALL GRIB_WRITE(IGRIB,IUNIT2) ! ! Free the memory. ! CALL GRIB_RELEASE(IGRIB) DEALLOCATE(VALUES_IN, VALUES_OUT) ! ! Loop back for next product. ! CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) END DO LOOP ! ! Close input and output files. ! CALL GRIB_CLOSE_FILE(IUNIT1, IRET) CALL GRIB_CLOSE_FILE(IUNIT2, IRET) ! STOP CONTAINS SUBROUTINE USAGE print*,'Usage: interpolation_example -i inputfile -o outputfile' print*, ' -i inputfile input file name (required)' print*, ' -o outputfile output file name (required)' RETURN END SUBROUTINE USAGE END PROGRAM INTF_EXAMPLE_WITH_ARRAY