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