1! 2! Copyright 1981-2016 ECMWF. 3! 4! This software is licensed under the terms of the Apache Licence 5! Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 6! 7! Unless required by applicable law or agreed to in writing, software 8! distributed under the License is distributed on an "AS IS" BASIS, 9! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 10! 11! In applying this licence, ECMWF does not waive the privileges and immunities 12! granted to it by virtue of its status as an intergovernmental organisation 13! nor does it submit to any jurisdiction. 14! 15PROGRAM EXAMPLE_INTF 16 17 USE grib_api 18 IMPLICIT NONE 19 INTEGER, DIMENSION(4) :: INTV 20 REAL, DIMENSION(4) :: REALV 21 CHARACTER(Len=20), DIMENSION(4) :: CHARV 22 23 INTEGER (KIND=4), DIMENSION(1) :: INGRIB(1), NEWFLD(1) 24 25 INTEGER :: IUNIT1, IUNIT2 26 INTEGER :: IGRIB 27 INTEGER :: IRET, NARGS 28 INTEGER (KIND=KINDOFINT) :: INLEN, NEWLEN 29 INTEGER :: J 30 INTEGER :: NRESOL 31 INTEGER :: Ni, Nj 32 INTEGER :: IBPV 33 INTEGER :: len_in, len_out 34 INTEGER :: IPARAM, ITABLE 35 REAL, DIMENSION(:), ALLOCATABLE :: VALUES_IN, VALUES_OUT 36 REAL :: DLAT, DLON 37 REAL :: VALMISS 38 39! Externals 40 INTEGER :: INTIN, INTOUT, INTF, IARGC 41 42 CHARACTER(Len=128) :: INFILE, OUTFILE 43 CHARACTER(Len=128), DIMENSION(20) :: CARG 44 CHARACTER(Len=128) :: MESSAGE = ' ' 45 CHARACTER(Len=32) :: CINGRIDT = ' ', CINGRIDN = ' ' 46 47! ********************************************************************** 48! Default values 49 50 IBPV = 16 ! bitsPerValue / accuracy 51 NRESOL = 0 ! Gaussian resolution 52 DLAT = 0.5 ! Latitude increment (in degrees) 53 DLON = 0.5 ! Longitude increment (in degrees) 54 VALMISS = -99999.9 ! Missing value indicator 55 56! Pick up file names from command line. 57 58 NARGS = 0 59 NARGS = IARGC() 60 IF ( NARGS .NE. 4) THEN 61 PRINT *, 'NARGS is incorrect' 62 CALL USAGE 63 STOP 64 END IF 65 66 DO J = 1, NARGS 67 CALL GETARG(J,CARG(J)) 68 END DO 69 70 DO J = 1, NARGS,2 71 IF (CARG(J) == '-i') THEN 72 INFILE=TRIM(CARG(J+1)) 73 ELSEIF (CARG(J) == '-o') THEN 74 OUTFILE=TRIM(CARG(J+1)) 75 ELSE 76 CALL USAGE 77 STOP 78 END IF 79 END DO 80 81 PRINT *, '*** Open input and output files' 82 83 CALL GRIB_OPEN_FILE (IUNIT1, INFILE, 'r', IRET) 84 IF ( IRET /= GRIB_SUCCESS) STOP ' GRIB_OPEN of INFILE failed' 85 86 CALL GRIB_OPEN_FILE (IUNIT2, OUTFILE, 'w', IRET) 87 IF ( IRET /= GRIB_SUCCESS) STOP ' GRIB_OPEN of OUTFILE failed' 88 89 PRINT *, '*** Start of loop through input GRIB-coded fields' 90 91 CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) 92 93 LOOP: DO WHILE (IRET /= GRIB_END_OF_FILE) 94 95 IF (IRET /= GRIB_SUCCESS) THEN 96 PRINT *, ' GRIB_NEW_FROM_FILE failed with error code = ', IRET 97 CALL GRIB_GET_ERROR_STRING(IRET,MESSAGE) 98 PRINT *, TRIM(MESSAGE) 99 PRINT *, 'Skipping to next message ...' 100 CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) 101 CYCLE 102 ENDIF 103 104 PRINT *, '*** Unpack the values in the input GRIB' 105 106 CALL GRIB_GET_SIZE(IGRIB,'values',len_in) 107 PRINT *, ' Allocating array for input field of size=',len_in 108 ALLOCATE(VALUES_IN(len_in)) 109 CALL GRIB_SET(IGRIB,'missingValue',VALMISS) 110 CALL GRIB_GET(IGRIB,'values',VALUES_IN) 111 112! Define the input field via calls to INTIN 113 CHARV(1) = 'unpacked' 114 IRET = INTIN('form',INTV, REALV, CHARV) 115 IF ( IRET /= 0 ) THEN 116 PRINT *, ' INTIN failed to set format' 117 STOP 118 ENDIF 119 120! Get parameter information (assumes unportable ParamId structure) 121 CALL GRIB_GET(IGRIB,'paramId',IPARAM) 122 ITABLE = 128 123 IF (IPARAM.GT.1000) THEN 124 ITABLE = IPARAM / 1000 125 IPARAM = MOD(IPARAM,1000) 126 ENDIF 127 128 INTV(1) = ITABLE 129 IRET = INTIN('table',INTV, REALV, CHARV) 130 IF ( IRET /= 0 ) THEN 131 PRINT *, ' INTIN failed to set table' 132 STOP 133 ENDIF 134 135 INTV(1) = IPARAM 136 IRET = INTIN('parameter',INTV, REALV, CHARV) 137 IF ( IRET /= 0 ) THEN 138 PRINT *, ' INTIN failed to set parameter' 139 STOP 140 ENDIF 141 142! Check it's a Gaussian grid and set resolution 143 CALL GRIB_GET(IGRIB,'gridType',CINGRIDT) 144 IF ((TRIM(CINGRIDT) == 'regular_gg') .OR. & 145 (TRIM(CINGRIDT) == 'reduced_gg')) THEN 146 CALL GRIB_GET(IGRIB,'gridName',CINGRIDN) 147 CHARV(1) = CINGRIDN(1:20) 148 PRINT *, ' gridName = ', CHARV(1) 149 IRET = INTIN('gridname',INTV, REALV, CHARV) 150 IF ( IRET /= 0 ) THEN 151 PRINT *, ' INTIN failed to set gridname' 152 STOP 153 ENDIF 154 ELSE 155 PRINT *, 'Input grid type ', TRIM(CINGRIDT), ' not regular_gg/reduced_gg' 156 CYCLE 157 ENDIF 158 159! Set the missing value 160 CHARV(1) = 'yes' 161 REALV(1) = VALMISS 162 IRET = INTIN('missingvalue',INTV, REALV, CHARV) 163 IF ( IRET /= 0 ) THEN 164 PRINT *, ' INTIN failed to set missingValue' 165 STOP 166 ENDIF 167 168! Define the output field with calls to INTOUT 169 CHARV(1) = 'unpacked' 170 IRET = INTOUT('form',INTV, REALV, CHARV) 171 IF ( IRET /= 0 ) THEN 172 PRINT *, ' INTOUT failed to set format' 173 STOP 174 ENDIF 175 176 PRINT *, '*** Define the grid interval for the new field(s)' 177 REALV(1) = DLAT 178 REALV(2) = DLON 179 IRET = INTOUT('grid', INTV, REALV, CHARV) 180 IF ( IRET /= 0 ) THEN 181 PRINT *, ' INTOUT failed to set grid resolution' 182 STOP 183 ENDIF 184 PRINT *, ' GRID TYPE = lat-lon' 185 PRINT *, ' RESOL = ', REALV(1), 'x', REALV(2) 186 187! Calculate the number of points (assumes global grid) 188 Ni = INT(360.0/REALV(2)) 189 Nj = INT(180.0/REALV(1)) + 1 190 len_out = Ni*Nj 191 PRINT *, ' Allocating output array of size = ', len_out 192 ALLOCATE(VALUES_OUT(LEN_OUT)) 193 194! Define the packing accuracy for the new field(s). 195 INTV(1) = IBPV 196 IRET = INTOUT('accuracy', INTV, REALV, CHARV) 197 IF ( IRET /= 0 ) THEN 198 PRINT *, ' INTOUT failed to set accuracy' 199 STOP 200 ENDIF 201 PRINT *, ' ACCURACY = ', IBPV 202 203 PRINT *, '*** Interpolate' 204 INLEN =1 205 NEWLEN = 1 206 IRET = INTF(INGRIB,len_in,VALUES_IN,NEWFLD,len_out,VALUES_OUT) 207 IF ( IRET /= 0 ) THEN 208 WRITE(*,*) ' INTF failed for product ' 209 WRITE(*,*) ' Skipping to next message ...' 210 CALL GRIB_RELEASE(IGRIB) 211 DEALLOCATE(VALUES_IN,VALUES_OUT) 212 CYCLE 213 ENDIF 214 215 PRINT *, '*** Set-up and write output GRIB' 216 217 CALL GRIB_SET(IGRIB,'gridType','regular_ll') 218 CALL GRIB_SET(IGRIB,'latitudeOfFirstGridPointInDegrees',90.0) 219 CALL GRIB_SET(IGRIB,'latitudeOfLastGridPointInDegrees',-90.0) 220 CALL GRIB_SET(IGRIB,'longitudeOfFirstGridPointInDegrees',0.0) 221 CALL GRIB_SET(IGRIB,'longitudeOfLastGridPointInDegrees',360.0 - DLON) 222 CALL GRIB_SET(IGRIB,'iDirectionIncrementInDegrees',DLON) 223 CALL GRIB_SET(IGRIB,'jDirectionIncrementInDegrees',DLON) 224 CALL GRIB_SET(IGRIB,'bitsPerValue',IBPV) 225 CALL GRIB_SET(IGRIB,'Ni',Ni) 226 CALL GRIB_SET(IGRIB,'Nj',Nj) 227 CALL GRIB_SET(IGRIB,'values',VALUES_OUT) 228 CALL GRIB_WRITE(IGRIB,IUNIT2) 229 230! Free the memory. 231 CALL GRIB_RELEASE(IGRIB) 232 DEALLOCATE(VALUES_IN, VALUES_OUT) 233 234! Loop back for next product. 235 CALL GRIB_NEW_FROM_FILE(IUNIT1, IGRIB, IRET) 236 237 END DO LOOP 238 239! Close input and output files. 240 CALL GRIB_CLOSE_FILE(IUNIT1, IRET) 241 CALL GRIB_CLOSE_FILE(IUNIT2, IRET) 242 243 STOP 244 245CONTAINS 246 247 SUBROUTINE USAGE 248 249 PRINT *, 'Usage: example_intf -i inputfile -o outputfile' 250 PRINT *, ' -i inputfile input file name (required)' 251 PRINT *, ' -o outputfile output file name (required)' 252 253 RETURN 254 255 END SUBROUTINE USAGE 256 257END PROGRAM EXAMPLE_INTF 258 259