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 OCEANP( INOCEAN, INLEN, OTOCEAN, OUTLEN) 12C 13C----> 14C**** OCEANP 15C 16C Purpose 17C ------- 18C 19C Interpolate GRIB format input ocean field to GRIB format ocean 20C field. 21C 22C 23C Interface 24C --------- 25C 26C IRET = OCEANP( INOCEAN, INLEN, OTOCEAN, OUTLEN) 27C 28C Input 29C ----- 30C 31C INOCEAN - Input ocean field (GRIB format). 32C INLEN - Input field length (words). 33C 34C 35C Output 36C ------ 37C 38C OTOCEAN - Output ocean field (GRIB format). 39C OUTLEN - Output field length (words). 40C 41C 42C Method 43C ------ 44C 45C Interpolate ocean field. 46C 47C If requested X grid increment does not match input GRIB 48C increment, do a full interpolation in both x and y. Otherwise, 49C interpolate only to make grid regular. 50C 51C Externals 52C --------- 53C 54C JDEBUG - Tests if debug output required. 55C GRIBEX - GRIB decoding/encoding. 56C INTLOG - Log error message. 57C JMEMHAN - Allocate/deallocate scratch memory. 58C INTOCN - Carry out ocean field to field interpolation. 59C 60C 61C Author 62C ------ 63C 64C J.D.Chambers ECMWF May 1996 65C 66C J.D.Chambers ECMWF Feb 1997 67C Allow for 64-bit pointers 68C 69C J.D.Chambers ECMWF Mar 1997 70C Allow for full or partial interpolation 71C (F or R option in call to INTOCN). 72C 73C----< 74C 75 IMPLICIT NONE 76C 77C Function arguments 78C 79 INTEGER INOCEAN, INLEN, OTOCEAN, OUTLEN 80C 81#include "parim.h" 82#include "nifld.common" 83#include "nofld.common" 84#include "intf.h" 85C 86C Parameters 87C 88 INTEGER JPROUTINE, JPALLOC, JPDEALL, JPSCR3, JPSCR4 89 INTEGER JP_LAT, JP_LONG, JP_GUESS 90 PARAMETER (JPROUTINE = 36900 ) 91 PARAMETER (JPALLOC = 1) 92 PARAMETER (JPDEALL = 0) 93 PARAMETER (JPSCR3 = 3) 94 PARAMETER (JPSCR4 = 4) 95 PARAMETER (JP_LONG = 3) 96 PARAMETER (JP_LAT = 4) 97 PARAMETER (JP_GUESS = 18049691) 98cs PARAMETER (JP_GUESS = 1675245) 99cs PARAMETER (JP_GUESS = 1237104) 100cs PARAMETER (JP_GUESS = 1038240) 101C `---> allow for 0.25*0.25 interpolated grid 102C 103C Local variables 104C 105 LOGICAL LHORIZN, LFLIP 106 CHARACTER*1 HFUNC 107 CHARACTER*1 HINTOPT 108 INTEGER IOLDN, IOLDW, IOLDS, IOLDE, IOLDEW, IOLDNS 109 INTEGER LOOPLAT, LOOPLON, I, J 110 REAL ZTEMP 111 INTEGER IERR, KPR, IWORD, IOLDSIZ, INEWSIZ, IINTPOL, ISTAGP 112 INTEGER NXP, NYP, INEWP, ITEMP 113 REAL RTOP, RBOTT, RRIGHT, RLEFT, RYRES, RXRES 114 REAL XARR, YARR 115 DIMENSION XARR(JPGRIB_ISEC1), YARR(JPGRIB_ISEC1) 116#ifndef _CRAYFTN 117#ifdef POINTER_64 118 INTEGER*8 INEW, IOLD 119#endif 120#endif 121 REAL NEW, OLD 122 POINTER ( INEW, NEW ) 123 POINTER ( IOLD, OLD ) 124 DIMENSION NEW( 1 ), OLD( 1 ) 125C 126C Externals 127C 128 INTEGER RESET_C, INTOCN, FIXAREA 129C 130C -----------------------------------------------------------------| 131C* Section 1. Initialise 132C -----------------------------------------------------------------| 133C 134 100 CONTINUE 135C 136C Test if debug output required. 137C 138 CALL JDEBUG() 139C 140 CALL INTLOG(JP_DEBUG, 141 X 'OCEANP: Trying to interpolate an ocean field.',JPQUIET) 142C 143 IF( .NOT.LNOGRID ) THEN 144 CALL INTLOG(JP_FATAL, 145 X 'OCEANP: GRID must be specified to interpolate an ocean field' 146 X , JPQUIET) 147 OCEANP = JPROUTINE + 1 148 GOTO 960 149 ENDIF 150C 151 OCEANP = 0 152 IERR = 0 153 KPR = 0 154 LFLIP = .FALSE. 155C 156C Save the original area and grid definitions 157C 158 IOLDN = NOAREA(1) 159 IOLDW = NOAREA(2) 160 IOLDS = NOAREA(3) 161 IOLDE = NOAREA(4) 162 IOLDEW = NOGRID(1) 163 IOLDNS = NOGRID(2) 164C 165C -----------------------------------------------------------------| 166C* Section 2. Unpack the input ocean field. 167C -----------------------------------------------------------------| 168C 169 200 CONTINUE 170C 171C Decode first four sections of GRIB code to find number of values 172C 173 IWORD = INLEN 174 IERR = 1 175 IOLDSIZ = OUTLEN 176 CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4, 177 X OLD, IOLDSIZ, INOCEAN, INLEN, IWORD, 'I',IERR) 178C 179 IF( IERR.NE.0 ) THEN 180 CALL INTLOG(JP_ERROR, 181 X 'OCEANP: GRIBEX option J decoding failed.',IERR) 182 OCEANP = IERR 183 GOTO 950 184 ENDIF 185C 186C Check it is an ocean field 187C 188 IF( (ISEC1(24).NE.1).OR.(ISEC1(37).NE.4) ) THEN 189 CALL INTLOG(JP_ERROR, 'OCEANP: Not an ocean field.',JPQUIET) 190 OCEANP = JPROUTINE + 3 191 GOTO 950 192 ENDIF 193C 194C Horizontal (lat/long) field GRIB setup different from others 195C 196 LHORIZN = (ISEC1(60).EQ.JP_LONG).AND.(ISEC1(61).EQ.JP_LAT) 197 IF( LHORIZN ) CALL INTLOG(JP_DEBUG, 198 X 'OCEANP: Horizontal (lat/long) field interpolation.',JPQUIET) 199C 200C Get scratch memory for unpacking the input ocean field. 201C 202 IOLDSIZ = ISEC2(2)*ISEC2(3) 203 CALL JMEMHAN( JPSCR4, IOLD, IOLDSIZ, JPALLOC, IERR) 204 IF( IERR.NE.0 ) THEN 205 CALL INTLOG(JP_ERROR, 206 X 'OCEANP: Scratch memory(4) allocation failed.',JPQUIET) 207 OCEANP = IERR 208 GOTO 950 209 ENDIF 210C 211C Decode data from GRIB code: -999.9 is a missing data value 212C 213 IWORD = INLEN 214 IERR = 1 215 ZSEC3(2) = -999.9 216 CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4, 217 X OLD, IOLDSIZ, INOCEAN, INLEN, IWORD, 'D',IERR) 218C 219 IF( (IERR.NE.0).AND.(IERR.NE.-2).AND.(IERR.NE.-4) ) THEN 220 CALL INTLOG(JP_ERROR,'OCEANP: GRIBEX decoding failed.',IERR) 221 OCEANP = IERR 222 GOTO 940 223 ENDIF 224C 225C -----------------------------------------------------------------| 226C* Section 3. Prepare the output ocean field. 227C -----------------------------------------------------------------| 228C 229 300 CONTINUE 230C 231C Setup interpolation options from input GRIB characteristics. 232C 233 IERR = RESET_C(ISEC1, ISEC2, ZSEC2, ISEC4) 234 IF( IERR.NE.0 ) THEN 235 CALL INTLOG(JP_ERROR, 236 X 'OCEANP: Setup interp. options from GRIB failed.',JPQUIET) 237 OCEANP = IERR 238 GOTO 940 239 ENDIF 240C 241C For horizontal field, match the desired area to the grid 242C 243 IERR = FIXAREA() 244 IF( IERR.NE.0 ) THEN 245 CALL INTLOG(JP_ERROR, 'OCEANP: FIXAREA failed.',JPQUIET) 246 OCEANP = IERR 247 GOTO 940 248 ENDIF 249C 250C Get scratch memory for output ocean field. 251C Have to guess how big the interpolated grid will be. 252C 253 INEWSIZ = JP_GUESS 254 CALL JMEMHAN(JPSCR3, INEW, INEWSIZ*2, JPALLOC, IERR) 255 IF( IERR.NE.0 ) THEN 256 CALL INTLOG(JP_ERROR, 257 X 'OCEANP: Scratch memory(3) allocation failed.',JPQUIET) 258 OCEANP = IERR 259 GOTO 940 260 ENDIF 261 INEWP = 1 262 ITEMP = 1 + INEWSIZ 263C 264C -----------------------------------------------------------------| 265C* Section 4. Interpolate ocean field. 266C -----------------------------------------------------------------| 267C 268 400 CONTINUE 269C 270 RTOP = FLOAT( NOAREA(1) ) / PPMULT 271 RBOTT = FLOAT( NOAREA(3) ) / PPMULT 272 RRIGHT = FLOAT( NOAREA(4) ) / PPMULT 273 RLEFT = FLOAT( NOAREA(2) ) / PPMULT 274C 275 RYRES = FLOAT( NOGRID(2) ) / PPMULT 276 RXRES = FLOAT( NOGRID(1) ) / PPMULT 277C 278C Select the interpolation option 279C (default = 'R' = interpolate only to make grid regular) 280C If requested X grid increment does not match input GRIB 281C increment, interpolate in x and y ('F' = full) 282C 283 HINTOPT = 'R' 284 IF( NOGRID(1).NE.ISEC2(9) ) HINTOPT = 'F' 285C 286 IERR = INTOCN(HINTOPT,'D', 287 X RTOP, RBOTT, RRIGHT, RLEFT, RYRES, RXRES, 288 X ISEC1, ISEC2, OLD, IOLDSIZ, 289 X INEWSIZ, NEW(ITEMP), XARR, YARR, NXP, NYP, 290 X NEW(INEWP), 291 X IINTPOL, ISTAGP) 292 IF( IERR.NE.0 ) THEN 293 CALL INTLOG(JP_ERROR, 'OCEANP: Interpolation failed.',JPQUIET) 294 OCEANP = IERR 295 GOTO 900 296 ENDIF 297C 298 NOAREA(1) = NINT( (YARR(1) + YARR(2)*(NYP-1)) * PPMULT ) 299 NOAREA(2) = NINT( XARR(1) * PPMULT ) 300 NOAREA(3) = NINT( YARR(1) * PPMULT ) 301 NOAREA(4) = NINT( (XARR(1) + XARR(2)*(NXP-1)) * PPMULT ) 302C 303 NOGRID(1) = NINT( ABS(XARR(2)) * PPMULT ) 304 NOGRID(2) = NINT( ABS(YARR(2)) * PPMULT ) 305 NONS = NYP 306 NOWE = NXP 307C 308C -----------------------------------------------------------------| 309C* Section 5. Code new field into GRIB. 310C -----------------------------------------------------------------| 311C 312 500 CONTINUE 313C 314C Re-order the rows from south to north to be north to south 315C (ECMWF convention). 316C 317 IF( LHORIZN.OR.(NOAREA(1).LT.NOAREA(3)) ) THEN 318 LFLIP = .TRUE. 319 DO LOOPLAT = 1, NYP/2 320 I = (LOOPLAT - 1)*NXP 321 J = (NYP - LOOPLAT)*NXP 322 DO LOOPLON = 1, NXP 323 I = I + 1 324 J = J + 1 325 ZTEMP = NEW(I) 326 NEW(I) = NEW(J) 327 NEW(J) = ZTEMP 328 ENDDO 329 ENDDO 330C 331 IF( .NOT.LHORIZN ) THEN 332 ITEMP = NOAREA(1) 333 NOAREA(1) = NOAREA(3) 334 NOAREA(3) = ITEMP 335C 336 ITEMP = ISEC2(4) 337 ISEC2(4) = ISEC2(7) 338 ISEC2(7) = ITEMP 339 ENDIF 340 ENDIF 341C 342C Setup GRIB sections for the output product: 343C 344C Use local definition 4 in GRIB section 1 ... 345C 346 ISEC1(73) = 0 347 ITEMP = 75 + ISEC1(71) + ISEC1(72) + ISEC1(73) + ISEC1(74) 348 ISEC1(ITEMP) = 0 349 350C 351 IF( LHORIZN ) THEN 352 ISEC1(62) = NOAREA(1)*10 353 ISEC1(63) = NOAREA(2)*10 354 ISEC1(64) = NOAREA(3)*10 355 ISEC1(65) = NOAREA(4)*10 356 ISEC1(66) = NOGRID(1)*10 357 ISEC1(67) = NOGRID(2)*10 358 IF( LFLIP ) ISEC1(67) = - ISEC1(67) 359 ISEC1(68) = 0 360 ISEC1(69) = 0 361C 362C Standard section 2 format for horizontal sections 363C 364 ISEC2(1) = 0 365 ISEC2(2) = NXP 366 ISEC2(3) = NYP 367 ISEC2(4) = NOAREA(1)/100 368 ISEC2(5) = NOAREA(2)/100 369 ISEC2(6) = 128 370 ISEC2(7) = NOAREA(3)/100 371 ISEC2(8) = NOAREA(4)/100 372 ISEC2(9) = NOGRID(1)/100 373 ISEC2(10) = NOGRID(2)/100 374Cjdc ISEC2(17) = 0 375 DO LOOPLAT = 11, 22 376 ISEC2(LOOPLAT) = 0 377 ENDDO 378 ELSE 379C Tim-Correct headers for vertical sections 380 IF( ISEC1(60).EQ.3 ) THEN 381 ISEC1(63) = NOAREA(2)*10 382 ELSEIF( ISEC1(60).EQ.4 ) THEN 383 ISEC1(63) = NOAREA(2)*10 384 ENDIF 385 IF( ISEC1(61).EQ.3 ) THEN 386 ISEC1(62) = NOAREA(1)*10 387 ELSEIF( ISEC1(61).EQ.4 ) THEN 388 ISEC1(62) = NOAREA(1)*10 389 ENDIF 390C 391 IF( ISEC1(61).EQ.1 ) THEN 392 ISEC1(64) = ISEC1(62) + NINT((NYP-1)*RYRES*1.0) 393 ISEC1(65) = ISEC1(63) + NINT((NXP-1)*RXRES*1000000.0) 394 ISEC1(66) = NOGRID(1)*NINT(1000000.0/PPMULT) 395 ISEC1(67) = NOGRID(2)/(PPMULT/1.0) 396 ELSE 397 ISEC1(64) = ISEC1(62) + NINT((NYP-1)*RYRES*1000.0) 398 ISEC1(65) = ISEC1(63) + NINT((NXP-1)*RXRES*1000000.0) 399 ISEC1(66) = NOGRID(1)*NINT(1000000.0/PPMULT) 400 ISEC1(67) = NOGRID(2)/(PPMULT/1000.0) 401 ENDIF 402C 403 IF( LFLIP ) ISEC1(67) = - ISEC1(67) 404 ISEC1(68) = 0 405C 406C Special ECMWF ocean short form of GRIB section 2 407C 408 ISEC2(1) = 192 409 ISEC2(2) = NXP 410 ISEC2(3) = NYP 411Cjdc ISEC2(4) = NOAREA(1)/100000 412Cjdc ISEC2(5) = NOAREA(2)/100 413Cjdc ISEC2(6) = 0 414Cjdc ISEC2(7) = NOAREA(3)/100000 415Cjdc ISEC2(8) = NOAREA(4)/100 416Cjdc DO LOOPLAT = 9, 22 417Cjdc ISEC2(LOOPLAT) = 0 418Cjdc ENDDO 419 ENDIF 420C 421C Ensure there is a section 3 in case there is missing data in the 422C output area 423C 424 ISEC1(5) = 192 425 ISEC3(1) = 0 426 ZSEC3(2) = -999.9 427C 428 ISEC4(1) = NXP * NYP 429C 430C If grid-point output, setup for 2nd order packing if requested. 431C 432 IF( (NOREPR.NE.JPSPHERE).AND.(NOREPR.NE.JPSPHROT) ) THEN 433 HFUNC = 'C' 434 IF( NOHFUNC.EQ.'K' ) THEN 435 HFUNC = 'K' 436 ISEC4(4) = 64 437 ISEC4(6) = 16 438 ISEC4(9) = 32 439 ISEC4(10) = 16 440 ISEC4(12) = 8 441 ISEC4(13) = 4 442 ISEC4(14) = 0 443 ISEC4(15) = -1 444 ELSE 445 ISEC4(4) = 0 446 ISEC4(6) = 0 447 ENDIF 448 ELSE 449 HFUNC = 'C' 450 IF( NOHFUNC.EQ.'C' ) THEN 451 ISEC2(6) = 2 452 ISEC4(4) = 64 453 ELSE IF( NOHFUNC.EQ.'S' ) THEN 454 ISEC2(6) = 1 455 ISEC4(4) = 0 456 ENDIF 457 ENDIF 458C 459C Switch off GRIBEX checking (sigh) 460C 461 CALL INTLOG(JP_DEBUG, 462 X 'OCEANP: switching OFF GRIBEX checking',JPQUIET) 463 CALL GRSVCK(0) 464C 465 IERR = 1 466 CALL GRIBEX( ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4, 467 X NEW,INEWSIZ,OTOCEAN,OUTLEN,IWORD,HFUNC,IERR) 468C 469 IF( IERR.NE.0 ) THEN 470 CALL INTLOG(JP_ERROR,'OCEANP: GRIBEX encoding failed.',IERR) 471 OCEANP = JPROUTINE + 5 472 GOTO 900 473 ENDIF 474 OUTLEN = IWORD 475C 476C Switch on GRIBEX checking (sigh again) 477C 478 CALL INTLOG(JP_DEBUG, 479 X 'OCEANP: switching ON GRIBEX checking',JPQUIET) 480 CALL GRSVCK(1) 481C 482C -----------------------------------------------------------------| 483C* Section 9. Closedown. 484C -----------------------------------------------------------------| 485C 486 900 CONTINUE 487C 488C Clear change flags for next product processing 489C 490 LCHANGE = .FALSE. 491 LSMCHNG = .FALSE. 492C 493C Return the scratch memory. 494C 495 CALL JMEMHAN( JPSCR3, INEW, INEWSIZ, JPDEALL, IERR) 496 IF ( IERR .NE. 0 ) THEN 497 CALL INTLOG(JP_ERROR, 498 X 'OCEANP: Scratch memory(3) reallocation failed.',JPQUIET) 499 OCEANP = IERR 500 ENDIF 501C 502 940 CONTINUE 503C 504 CALL JMEMHAN( JPSCR4, IOLD, IOLDSIZ, JPDEALL, IERR) 505 IF ( IERR .NE. 0 ) THEN 506 CALL INTLOG(JP_ERROR, 507 X 'OCEANP: Scratch memory(4) reallocation failed.',JPQUIET) 508 OCEANP = IERR 509 ENDIF 510C 511 950 CONTINUE 512C 513C Restore the original area and grid definitions 514C 515 NOAREA(1) = IOLDN 516 NOAREA(2) = IOLDW 517 NOAREA(3) = IOLDS 518 NOAREA(4) = IOLDE 519 NOGRID(1) = IOLDEW 520 NOGRID(2) = IOLDNS 521C 522 960 CONTINUE 523C 524 CALL INTLOG(JP_DEBUG, 525 X 'OCEANP: Returning from interpolating an ocean field.',JPQUIET) 526C 527 RETURN 528 END 529