1C Copyright 1981-2016 ECMWF.
2C
3C This software is licensed under the terms of the Apache Licence
4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5C
6C In applying this licence, ECMWF does not waive the privileges and immunities
7C granted to it by virtue of its status as an intergovernmental organisation
8C nor does it submit to any jurisdiction.
9C
10
11       INTEGER FUNCTION INTF( INGRIB,INLEN,FLDIN,OUTGRIB,OUTLEN,FLDOUT)
12C
13C---->
14C**** INTF
15C
16C     Purpose
17C     -------
18C
19C     Interpolate input field...
20C
21C
22C     Interface
23C     ---------
24C
25C     IRET = INTF( INGRIB,INLEN,FLDIN,OUTGRIB,OUTLEN,FLDOUT)
26C
27C     Input
28C     -----
29C
30C     INGRIB - Input field (packed).
31C     INLEN  - Input field length (words).
32C     FLDIN  - Input field (unpacked).
33C
34C
35C     Output
36C     ------
37C
38C     OUTGRIB - Output field (packed).
39C     OUTLEN  - Output field length (words).
40C     FLDOUT  - Output field (unpacked).
41C
42C
43C     Method
44C     ------
45C
46C     Call interpolation routines.
47C
48C
49C     Externals
50C     ---------
51C
52C     IBASINI - Ensure basic interpolation setup is done.
53C     INSANE  - Ensure no outrageous values given for interpolation.
54C     PDDEFS  - Setup interpolation using parameter dependent options.
55C     PRECIP  - Says if field is to have 'precipitation' treatment
56C     INTFAU  - Prepare to interpolate unpacked input field.
57C     INTFAP  - Prepare to interpolate packed input field.
58C     HNTFAU  - Prepare to interpolate unpacked input field (using
59C               Hirlam 12-point for rotations).
60C     HNTFAP  - Prepare to interpolate packed input field (using
61C               Hirlam 12-point for rotations).
62C     INTFB   - Interpolate input field.
63C     INTLOG  - Log error message.
64C     JDEBUG  - Checks environment to switch on/off debug
65C     INTWAVE - Interpolate quasi-regular lat/long wave field
66C               to a regular lat/long field.
67C     OCEANP  - Interpolate GRIB ocean field.
68C     RESET_C - Reset interpolation handling options using GRIB product.
69C
70C
71C     Author
72C     ------
73C
74C     J.D.Chambers     ECMWF     Aug 1994
75C
76C----<
77C
78      IMPLICIT NONE
79C
80C     Function arguments
81C
82      INTEGER INGRIB(*),OUTGRIB(*),INLEN,OUTLEN
83      REAL FLDIN(*),FLDOUT(*)
84C
85#include "parim.h"
86#include "nifld.common"
87#include "nofld.common"
88#include "grfixed.h"
89#include "intf.h"
90#include "current.h"
91C
92C     Parameters
93C
94      INTEGER JPROUTINE
95      PARAMETER (JPROUTINE = 26000 )
96C
97C     Local variables
98C
99      INTEGER IWORD, IRET, IORIGLN, ISAME
100      INTEGER LOOP, IIHOLD, IOHOLD
101      DIMENSION IIHOLD(4), IOHOLD(4)
102      INTEGER NUMTABL, NUMPROD
103      LOGICAL L98WAVE, LUNROT
104      CHARACTER*1 HTYPE
105C
106C     Externals
107C
108      INTEGER IBASINI, PDDEFS, INSANE
109      INTEGER INTFAU, INTFAP, INTFB, HNTFAU, HNTFAP
110      INTEGER INTWAVE2, INTWAVU, OCEANP, OCEANU, RESET_C !, INTWAVE
111      LOGICAL PRECIP
112C
113C     -----------------------------------------------------------------|
114C*    Section 1.   Initialise
115C     -----------------------------------------------------------------|
116C
117  100 CONTINUE
118C
119      INTF = 0
120      IRET = 0
121      IORIGLN = OUTLEN
122      NOMISS = 0
123      LUNROT = .FALSE.
124      OUTLROT = 0
125
126C
127C     Save input and output area definitions
128C
129      DO 110 LOOP = 1, 4
130        IIHOLD(LOOP) = NIAREA(LOOP)
131        IOHOLD(LOOP) = NOAREA(LOOP)
132  110 CONTINUE
133C
134C     Check if debug option turned on
135C
136      CALL JDEBUG()
137C
138C     If the input is a set of U and V rotated lat/long fields,
139C     set return length to zero and do not copy input to output
140C
141      IF( (NIREPR.EQ.JPREGULAR).AND.LWIND.AND.LNOROTA ) THEN
142        OUTLEN = 0
143        IRET   = 0
144        CALL INTLOG(JP_DEBUG,
145     X    'INTF: Input U and V rotated lat/long fields ...',JPQUIET)
146        CALL INTLOG(JP_DEBUG,
147     X    'INTF: ... no further interpolation has been done',JPQUIET)
148        GOTO 900
149      ENDIF
150C
151C     -----------------------------------------------------------------|
152C*    Section 2.   Prepare to interpolate input field.
153C     -----------------------------------------------------------------|
154C
155  200 CONTINUE
156C
157C     Ensure that basic initialisation has been done
158C
159      IRET = IBASINI(0)
160      IF( IRET.NE.0 ) THEN
161        CALL INTLOG(JP_ERROR,'INTF: basic initialisation fail.',JPQUIET)
162        INTF = IRET
163        GOTO 900
164      ENDIF
165C
166C     Allocate work array ZNFELDI if not already done.
167C
168       LACCUR = .FALSE.
169      IF( NIFORM.EQ.1 ) THEN
170        IF( IZNJDCI.NE.1952999238 ) THEN
171          CALL JMEMHAN( 19, IZNFLDI, JPEXPAND, 1, IRET)
172          IF( IRET.NE.0 ) THEN
173            CALL INTLOG(JP_WARN,'INTF: ZNFELDI allocation fail',JPQUIET)
174            INTF = IRET
175            GOTO 900
176          ENDIF
177          IZNJDCI = 1952999238
178        ENDIF
179C
180C       Unpack the field headers for packed input.
181C
182        CALL INTLOG(JP_DEBUG, 'INTF: Unpack GRIB headers.',JPQUIET)
183        IRET = 1
184        IWORD = 0
185        CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4,
186     X              ZNFELDI, JPEXPAND, INGRIB, INLEN, IWORD, 'J',IRET)
187C
188        IF( (IRET.NE.0).AND.(IRET.NE.811) ) THEN
189          CALL INTLOG(JP_ERROR,
190     X      'INTF: Failed to unpack GRIB heders.',JPQUIET)
191          INTF = IRET
192          GOTO 900
193        ENDIF
194C
195C       Reset interpolation handling options using GRIB values.
196C
197        IRET = RESET_C( ISEC1, ISEC2, ZSEC2, ISEC4)
198        IF( IRET.NE.0 ) THEN
199          CALL INTLOG(JP_WARN,
200     X      'INTF: Setup of interp. options from GRIB failed',JPQUIET)
201          INTF = IRET
202          GOTO 900
203        ENDIF
204C end of NIFORM = 1
205       LACCUR = .TRUE.
206      ENDIF
207C
208CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
209cs   Regular Gaussian has to be set here for unpacked fields
210      IF( NIFORM.NE.1 ) THEN
211cs    this is for merging with grib_api
212        LUVCOMP = .FALSE.
213        IF( NOREPR.EQ.JPNOTYPE ) THEN
214         IF( (NOGAUSO.NE.NOGAUSS).OR.(HOGAUST.NE.'F') ) THEN
215            HTYPE = 'F'
216            CALL JGETGG(NOGAUSS,HTYPE,ROGAUSS,NOLPTS,IRET)
217            IF( IRET.NE.0 ) THEN
218              CALL INTLOG(JP_ERROR,
219     X          'INTF: JGETGG failed, NOGAUSS = ',NOGAUSS)
220              INTF = IRET
221              GOTO 900
222            ENDIF
223            NOGAUSO = NOGAUSS
224            HOGAUST = HTYPE
225          ENDIF
226          NOREPR = JPGAUSSIAN
227        ENDIF
228      ENDIF
229CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
230
231C     Check that no outrageous values given for interpolation
232C
233      ISAME = INSANE()
234      IF ( ISAME .GT. 0 ) THEN
235        CALL INTLOG(JP_ERROR,
236     X    'INTF: Interpolation cannot use given values.',JPQUIET)
237        INTF = ISAME
238        GOTO 900
239      ENDIF
240Cs setting of out Date
241      NODATE = NIDATE
242C
243C     If output is same as the input, set return length to zero and
244C     do not copy input to output
245C
246      IF( ISAME.EQ.-1 ) THEN
247        OUTLEN = 0
248        IRET   = 0
249        CALL INTLOG(JP_DEBUG,
250     X    'INTF: Output is same as the input.',JPQUIET)
251        CALL INTLOG(JP_DEBUG,
252     X    'INTF: No interpolation carried out.',JPQUIET)
253        GOTO 900
254      ENDIF
255C
256C     Set precipitation flag if user hasn't
257C
258      IF( .NOT.LPRECSET ) LPREC = PRECIP()
259C
260C     Handle packed fields
261C
262C
263C     -----------------------------------------------------------------|
264C*    Section 3.   Use special interpolation for:
265C                  - an ECMWF wave field.
266C                  - a reduced latitude-longitude field
267C     -----------------------------------------------------------------|
268C
269  300   CONTINUE
270C
271        NUMTABL = NITABLE
272        NUMPROD = NUMTABL*1000 + NIPARAM
273        L98WAVE = (NUMTABL.EQ.140).OR.
274     X            (NUMPROD.EQ.131229).OR.
275     X            (NUMPROD.EQ.131232).OR.
276     X            (NIREPR.EQ.26)
277      IF( L98WAVE ) THEN
278          CALL INTLOG(JP_DEBUG,
279     X      'INTF: Wave-type interpolation required.',JPQUIET)
280C
281        OUTLEN = IORIGLN
282        IF( NIFORM.EQ.1 ) THEN
283             IRET = INTWAVE2(INGRIB,INLEN,OUTGRIB,OUTLEN)
284cs             IRET = INTWAVE(INGRIB,INLEN,OUTGRIB,OUTLEN)
285        ELSE
286            IRET = INTWAVU(FLDIN,INLEN,FLDOUT,OUTLEN)
287        ENDIF
288          IF( IRET.EQ.0 ) THEN
289            CALL INTLOG(JP_DEBUG,
290     X        'INTF: Wave-type interpolated OK.',JPQUIET)
291          ELSE
292            CALL INTLOG(JP_DEBUG,
293     X        'INTF: Wave-type interpolation failed.',JPQUIET)
294          ENDIF
295          INTF = IRET
296          GOTO 900
297        ENDIF
298C
299C     -----------------------------------------------------------------|
300C*    Section 4.   Use special interpolation for an ECMWF ocean field.
301C     -----------------------------------------------------------------|
302C
303  400   CONTINUE
304C
305        IF( ((ISEC1(24).EQ.1).AND.(ISEC1(37).EQ.4)).OR.LOCEAN ) THEN
306C
307          CALL INTLOG(JP_DEBUG,
308     X      'INTF: Ocean field interpolation required.',JPQUIET)
309C
310          OUTLEN = IORIGLN
311          IF( NIFORM.EQ.1 ) THEN
312            IRET = OCEANP(INGRIB,INLEN,OUTGRIB,OUTLEN)
313          ELSE
314            IRET = OCEANU(FLDIN,INLEN,FLDOUT,OUTLEN)
315          ENDIF
316          IF( IRET.EQ.0 ) THEN
317            CALL INTLOG(JP_DEBUG,
318     X        'INTF: Ocean field interpolated OK.',JPQUIET)
319          ELSE
320            CALL INTLOG(JP_DEBUG,
321     X        'INTF: Ocean field interpolation failed.',JPQUIET)
322          ENDIF
323          INTF = IRET
324          GOTO 900
325        ENDIF
326C
327C     -----------------------------------------------------------------|
328C*    Section 7.   Continue interpolation setup in other cases.
329C     -----------------------------------------------------------------|
330C
331  500   CONTINUE
332      IF( NIFORM.EQ.1 ) THEN
333C
334C       Unpack (and rotate) field if necessary
335C
336        IF( LUSEHIR ) THEN
337          IRET = HNTFAP(INGRIB,INLEN)
338        ELSE
339          IRET = INTFAP(INGRIB,INLEN)
340        ENDIF
341C
342C       If a bitmap encountered with some missing values (IRET=-4),
343C       product cannot be interpolated unless 'missingvalue'
344C       specified via INTIN
345C
346        IF( .NOT.LIMISSV ) THEN
347          IF( (IRET.NE.0).AND.(IRET.NE.-2) ) THEN
348            IF( IRET.EQ.-4 ) THEN
349              CALL INTLOG(JP_WARN,
350     X        'INTF: Product has bitmap and missing data.',JPQUIET)
351              CALL INTLOG(JP_WARN,
352     X        'INTF: Try Using INTIN "missingvalue" option',JPQUIET)
353            ENDIF
354            INTF = -4
355            GOTO 900
356          ENDIF
357        ELSE
358          IF( IRET.GT.0 ) THEN
359            CALL INTLOG(JP_WARN,
360     X        'INTF: Problems preparing for interpolation.',JPQUIET)
361            INTF = IRET
362            GOTO 900
363          ENDIF
364        ENDIF
365C
366C     Handle unpacked fields
367C
368      ELSE
369        LUNROT = .TRUE.
370        IF( LUSEHIR ) THEN
371          IRET = HNTFAU(FLDIN,INLEN)
372        ELSE
373          IRET = INTFAU(FLDIN,INLEN)
374        ENDIF
375      ENDIF
376C
377      IF( (IRET.NE.0).AND.(IRET.NE.-2) ) THEN
378        CALL INTLOG(JP_WARN,'INTF: Prepare interpolate fail',JPQUIET)
379        INTF = IRET
380        GOTO 900
381      ENDIF
382C
383C     Field values are now in ZNFELDI.
384C
385C     Setup output length same as input GRIB length in case straight
386C     copy is done later (ie input is transferred direct to output
387C     without postprocessing).
388C
389      OUTLEN = INLEN
390
391C
392C     Setup interpolation options based on parameter in field.
393C
394      IRET = PDDEFS()
395      IF( IRET.NE.0 ) THEN
396        CALL INTLOG(JP_ERROR,
397     X    'INTF: Setup interpolation options from param failed',JPQUIET)
398        INTF = IRET
399        GOTO 900
400      ENDIF
401C
402C     -----------------------------------------------------------------|
403C*    Section 8.   Interpolate input field.
404C     -----------------------------------------------------------------|
405C
406  700 CONTINUE
407C
408C     If all values missing, set flag to ensure all are missing in the
409C     interpolated field.
410C
411      IF( ISEC4(1).LT.0 ) NOMISS = ISEC4(1)
412C
413C     Perform the interpolation.
414C
415      CALL INTLOG(JP_DEBUG,'INTF: Perform the interpolation.',JPQUIET)
416C
417      OUTLEN = IORIGLN
418      IF( NIFORM.EQ.1 ) THEN
419        IRET = INTFB( INGRIB,INLEN,OUTGRIB,OUTLEN,FLDOUT)
420      ELSE
421        IRET = INTFB( ZNFELDI,INLEN,OUTGRIB,OUTLEN,FLDOUT)
422      ENDIF
423C
424      IF( LUNROT.AND.LNOROTA ) THEN
425          OUTLEN = OUTLROT
426      ENDIF
427C
428      IF( IRET.NE.0 ) THEN
429        CALL INTLOG(JP_ERROR,'INTF: Interpolation failed.',JPQUIET)
430        INTF = IRET
431        GOTO 900
432      ELSE
433C
434        CALL INTLOG(JP_DEBUG,
435     X    'INTF: Interpolation finished successfully.',JPQUIET)
436      ENDIF
437C
438C     -----------------------------------------------------------------|
439C*    Section 9.   Closedown.
440C     -----------------------------------------------------------------|
441C
442  900 CONTINUE
443C
444C     Clear change flags for next product processing
445C
446      LCHANGE = .FALSE.
447      LSMCHNG = .FALSE.
448C
449C     Restore input and output area definitions
450C
451      DO 910 LOOP = 1, 4
452        NIAREA(LOOP) = IIHOLD(LOOP)
453        NOAAPI(LOOP) = NOAREA(LOOP)
454        NOAREA(LOOP) = IOHOLD(LOOP)
455  910 CONTINUE
456C
457      RETURN
458      END
459