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