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 INTWAVE2(INGRIB,INLEN,OUTGRIB,OUTLEN) 12C 13C----> 14C**** INTWAVE2 15C 16C Purpose 17C ------- 18C 19C Interpolate an ECMWF wave field wave field or a reduced 20C latitude-longitude field 21C 22C 23C Interface 24C --------- 25C 26C IRET = INTWAVE2(INGRIB,INLEN,OUTGRIB,OUTLEN) 27C 28C Input 29C ----- 30C 31C INGRIB - Input field (packed). 32C INLEN - Input field length (words). 33C 34C 35C Output 36C ------ 37C 38C OUTGRIB - Output field (packed). 39C OUTLEN - Output field length (words). 40C 41C 42C Method 43C ------ 44C 45C Call interpolation routine; then repack into GRIB format. 46C 47C 48C Externals 49C --------- 50C 51C FIXAREA - Fixup area definition to correspond to grid definitions 52C GRIBEX - Decode/encode GRIB product. 53C JMEMHAN - Handles memory allocation. 54C JDEBUG - Checks environment to switch on/off debug 55C INTLOG - Log error message. 56C WAVEXX1 - Interpolate wave fields (except 2D spectra) 57C WV2DXXX - Interpolate wave 2D spectra fields 58C MKFRAME - Create a 'frame' from a rectangular field 59C MKBITMP - Apply a bitmap to a rectangular field 60C 61C 62C Author 63C ------ 64C 65C S.Curic ECMWF Jun 2009 66C 67C 68C----< 69C 70 IMPLICIT NONE 71C 72C Function arguments 73 INTEGER INLEN,OUTLEN 74 INTEGER INGRIB(*),OUTGRIB(*) 75C 76#include "parim.h" 77#include "nifld.common" 78#include "nofld.common" 79#include "intf.h" 80#include "grfixed.h" 81#include "jparams.h" 82C 83C Parameters 84 INTEGER JPROUTINE, JP2DSP, JP2DSPQ, JPMAXLT 85 86 PARAMETER (JPROUTINE = 40200 ) 87 PARAMETER (JP2DSP = 250 ) ! for 2D wave spectra (whole) 88 PARAMETER (JP2DSPQ = 251 ) ! for 2D wave spectra (single) 89 PARAMETER (JPMAXLT = 1801 ) 90C 91C Local variables 92C 93 CHARACTER*1 HFUNC 94 95 REAL ZMISS 96 DATA ZMISS/-9999999.0/ 97 INTEGER IWORD, IRET, ISIZE, NSPEC, LAT1, LOOP, IERR 98 INTEGER KNUM, NUM_E_W, NUM_N_S, NFULLNS, NFULLEW 99 REAL GRIDWE, GRIDNS, NORTH, SOUTH, WEST, EAST 100 INTEGER NUMPTS 101 DIMENSION NUMPTS(JPMAXLT) 102#ifndef _CRAYFTN 103#ifdef POINTER_64 104 INTEGER*8 INEWAVE, IDISTNW, IZ2DSP, INEWIDX 105#endif 106#endif 107 REAL NEWAVE 108 POINTER ( INEWAVE, NEWAVE ) 109 DIMENSION NEWAVE( 1 ) 110 REAL*4 DISTNEW 111 POINTER ( IDISTNW, DISTNEW ) 112 DIMENSION DISTNEW( 1 ) 113 REAL Z2DSP 114 POINTER ( IZ2DSP, Z2DSP ) 115 DIMENSION Z2DSP( 1 ) 116 INTEGER NEWIDX 117 POINTER ( INEWIDX, NEWIDX ) 118 DIMENSION NEWIDX( 1 ) 119 REAL RNLAT, RSLAT 120 INTEGER NUMTABL, NUMPROD 121 LOGICAL LCOEFFS, L98WAVE, LFRAME, LBITMP 122 LOGICAL LDEBUG 123 INTEGER NPARAM 124 INTEGER MISSLAT, KOLDNUM 125C 126C Externals 127 INTEGER WAVEXX2, WV2DXX2, AREACHK, FIXAREA 128 INTEGER NUMPTNS, NUMPTWE, JNORSGG, MKBITMP 129C 130C ------------------------------------------------------------------ 131C Section 1. Initialise 132C ------------------------------------------------------------------ 133C 134 100 CONTINUE 135C 136 INTWAVE2 = 0 137C 138 LCOEFFS = .FALSE. 139C 140C Check if debug option turned on 141 CALL JDEBUG() 142 LDEBUG = NDBG.GT.0 143C 144C Allocate work array ZNFELDI if not already done. 145C 146 IF( IZNJDCI.NE.1952999238 ) THEN 147 CALL JMEMHAN( 19, IZNFLDI, JPEXPAND, 1, IRET) 148 IF( IRET.NE.0 ) THEN 149 CALL INTLOG(JP_WARN,'INTWAVE2: ZNFELDI allocate fail',JPQUIET) 150 INTWAVE2= IRET 151 GOTO 900 152 ENDIF 153 IZNJDCI = 1952999238 154 ENDIF 155C 156C Unpack GRIB message headers (using ZMISS as missing data value). 157C 158cs IWORD = INLEN 159 IRET = 1 160 ISEC3(2) = NINT(ZMISS) 161 ZSEC3(2) = ZMISS 162 CALL GRIBEX (ISEC0,ISEC1,ISEC2,ZSEC2,ISEC3,ZSEC3,ISEC4, 163 X ZNFELDI,JPEXPAND,INGRIB,INLEN,IWORD,'D',IRET) 164 165 IF(IRET.GT.0) THEN 166 CALL INTLOG(JP_ERROR,'INTWAVE2: GRIBEX decoding failed.',JPQUIET) 167 INTWAVE2 = JPROUTINE + 1 168 GOTO 900 169 ENDIF 170 171 NPARAM = ISEC1(6) 172C 173C Check that the field can be handled 174C 175 NUMTABL = ISEC1(2)*1000 + ISEC1(1) 176 NUMPROD = NUMTABL*1000 + ISEC1(6) 177 L98WAVE = (NUMTABL.EQ.98140).OR. 178 X (NUMPROD.EQ.98131229).OR. 179 X (NUMPROD.EQ.98131232).OR. 180 X (NIREPR.EQ.26) 181 IF( .NOT.L98WAVE ) THEN 182 CALL INTLOG(JP_WARN,'INTWAVE2: Not an ECMWF wave field',JPQUIET) 183 INTWAVE2 = JPROUTINE + 2 184 GOTO 900 185 ENDIF 186 IF( ISEC2(1).NE.0 ) THEN 187 CALL INTLOG(JP_WARN, 188 X 'INTWAVE2: Not a lat/long field',JPQUIET) 189 INTWAVE2 = JPROUTINE + 3 190 GOTO 900 191 ENDIF 192 IF( (ISEC2(6).NE.128) .OR. (ISEC2(10).LT.10) ) THEN 193 CALL INTLOG(JP_ERROR, 194 X 'INTWAVE2: Cannot handle longitude increment',ISEC2(10)) 195 INTWAVE2 = JPROUTINE + 4 196 GOTO 900 197 ENDIF 198C 199C ------------------------------------------------------------------ 200C Section 2. Sort out interpolation area. 201C ------------------------------------------------------------------ 202C 203 200 CONTINUE 204C 205 NSPEC = 1 206 KNUM = ISEC2(3) 207C 208C Adjust (sub-)area limits to suit the grid 209C 210 IRET = FIXAREA() 211 IF ( IRET .NE. 0 ) THEN 212 CALL INTLOG(JP_ERROR, 'INTWAVE2: FIXAREA failed.',JPQUIET) 213 INTWAVE2 = JPROUTINE + 5 214 GOTO 900 215 ENDIF 216 NORTH = FLOAT(NOAREA(1))/PPMULT 217 WEST = FLOAT(NOAREA(2))/PPMULT 218 SOUTH = FLOAT(NOAREA(3))/PPMULT 219 EAST = FLOAT(NOAREA(4))/PPMULT 220C 221C Calculate the number of points E-W and N-S in new grid area 222C 223 IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN)) THEN 224 GRIDWE = NOGAUSS 225 GRIDNS = 0 226 ELSE IF( NOREPR.EQ.JPREDLL ) THEN 227 GRIDWE = 360.0/NOLPTS((NOREDLL/2)+1) 228 GRIDNS = FLOAT(NOGRID(2))/PPMULT 229 ELSE 230 GRIDWE = FLOAT(NOGRID(1))/PPMULT 231 GRIDNS = FLOAT(NOGRID(2))/PPMULT 232 ENDIF 233C 234 IRET = AREACHK(GRIDWE,GRIDNS,NORTH,WEST,SOUTH,EAST) 235 IF( IRET.NE.0 ) THEN 236 INTWAVE2 = JPROUTINE + 6 237 GOTO 900 238 ENDIF 239 NOAREA(1) = NINT(NORTH*PPMULT) 240 NOAREA(2) = NINT(WEST*PPMULT) 241 NOAREA(3) = NINT(SOUTH*PPMULT) 242 NOAREA(4) = NINT(EAST*PPMULT) 243C 244 IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN)) THEN 245 GRIDWE = 360.0/NOLPTS(NOGAUSS) 246 NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE) 247 NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1 248 NUM_N_S = JNORSGG(SOUTH,ROGAUSS,NOGAUSS,1) - 249 X JNORSGG(NORTH,ROGAUSS,NOGAUSS,0) + 1 250 NFULLNS = NOGAUSS*2 251 ELSE IF( NOREPR.EQ.JPREDLL ) THEN 252 NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE) 253 NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1 254 NUM_N_S = NUMPTNS(NORTH,SOUTH,GRIDNS) 255 NFULLNS = NOREDLL 256 ELSE 257 NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE) 258 NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1 259 NUM_N_S = NUMPTNS(NORTH,SOUTH,GRIDNS) 260 NFULLNS = NUM_N_S 261 ENDIF 262C 263C 264C-------------------------------------------------------------------- 265 KNUM = NINS 266cs KNUM = ISEC2(3) 267 IF( ISEC2(17).EQ.1 ) THEN 268C 269C Input field is a reduced latitude/longitude grid 270C 271C .. but it may be 'pseudo-gaussian' in layout 272C (ie global, symmetric about the equator but no latitude 273C at the equator) 274C 275 IF( (ISEC2(4).NE.90000).AND.(MOD(ISEC2(3),2).EQ.0) ) THEN 276C 277 DO LOOP = 1, ISEC2(3) 278 NUMPTS(LOOP) = ISEC2(22+LOOP) 279 ENDDO 280C 281 ELSE 282C 283 MISSLAT = (90000 - ISEC2(4))/ISEC2(10) 284 DO LOOP = 1, MISSLAT 285 NUMPTS(LOOP) = 0 286 ENDDO 287 KOLDNUM = 1 + (90000 - ISEC2(7))/ISEC2(10) 288 DO LOOP = 1, (KOLDNUM-MISSLAT) 289 NUMPTS(LOOP+MISSLAT) = ISEC2(22+LOOP) 290 ENDDO 291 DO LOOP = (KOLDNUM+1), KNUM 292 NUMPTS(LOOP) = 0 293 ENDDO 294 ENDIF 295 ENDIF 296C 297C 298C ------------------------------------------------------------------ 299C Section 3. Interpolate wave fields other than 2D spectra. 300C ------------------------------------------------------------------ 301C 302 300 CONTINUE 303C 304C Handle if not 2D spectra .. 305C 306 IF( (ISEC1(6).NE.JP2DSP).AND.(ISEC1(6).NE.JP2DSPQ) ) THEN 307 IF( LDEBUG ) THEN 308 CALL INTLOG(JP_DEBUG, 309 X 'INTWAVE2: Interpolate wave flds other than 2D spectra',JPQUIET) 310 ENDIF 311C 312C Get some scratch space for the interpolated field. 313C 314 ISIZE = NFULLEW * NUM_N_S 315 CALL JMEMHAN( 3, INEWAVE, ISIZE, 1, IRET) 316 IF ( IRET.NE.0 ) THEN 317 CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET) 318 INTWAVE2 = JPROUTINE + 7 319 GOTO 900 320 ENDIF 321C 322 ISIZE = NFULLEW * NFULLNS * 4 323 CALL JMEMHAN( 4, INEWIDX, ISIZE, 1, IRET) 324 IF ( IRET.NE.0 ) THEN 325 CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET) 326 INTWAVE2 = JPROUTINE + 8 327 GOTO 900 328 ENDIF 329C 330 ISIZE = NFULLEW * NFULLNS * 10 331 CALL JMEMHAN( 5, IDISTNW, ISIZE, 1, IRET) 332 IF ( IRET.NE.0 ) THEN 333 CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET) 334 INTWAVE2 = JPROUTINE + 9 335 GOTO 900 336 ENDIF 337 338C Interpolate the field 339C 340 IRET = WAVEXX2(NPARAM,KNUM,NUMPTS,NUM_N_S,GRIDNS,GRIDWE,ZNFELDI, 341 X NEWAVE,NORTH,WEST,ZMISS) 342 IF( IRET.NE.0 ) THEN 343 CALL INTLOG(JP_ERROR,'INTWAVE2: Interpolation failed.',JPQUIET) 344 INTWAVE2 = JPROUTINE + 10 345 GOTO 900 346 ENDIF 347C 348C ------------------------------------------------------------------ 349C Section 4. Interpolate wave 2D spectra field. 350C ------------------------------------------------------------------ 351C 352 ELSE 353C 354 400 CONTINUE 355 356 IF( LDEBUG ) THEN 357 CALL INTLOG(JP_DEBUG, 358 X 'INTWAVE2: Interpolate wave 2D spectra field',JPQUIET) 359 ENDIF 360C 361C Find number of 2D spectra values at each point 362C 363 IF( ISEC4(8).NE.0 ) LCOEFFS = .TRUE. 364 IF( (ISEC4(8).NE.0).AND.(ISEC1(6).NE.JP2DSPQ) ) THEN 365 NSPEC = ISEC4(50) * ISEC4(51) 366 NIMATR = 1 367 ELSE 368 NSPEC = 1 369 NIMATR = 0 370 ENDIF 371C 372C Get some scratch space for the interpolated field. 373C 374 ISIZE = NSPEC * NFULLEW * NUM_N_S 375 CALL JMEMHAN( 3, INEWAVE, ISIZE, 1, IRET) 376 IF ( IRET.NE.0 ) THEN 377 CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET) 378 INTWAVE2 = JPROUTINE + 11 379 GOTO 900 380 ENDIF 381C 382 ISIZE = NFULLEW * NFULLNS 383 CALL JMEMHAN( 4, INEWIDX, ISIZE, 1, IRET) 384 IF ( IRET.NE.0 ) THEN 385 CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET) 386 INTWAVE2 = JPROUTINE + 12 387 GOTO 900 388 ENDIF 389C 390C Interpolate the field 391C 392 IRET = WV2DXX2(NIMATR,KNUM,NUMPTS,NUM_N_S,GRIDNS,GRIDWE,ZNFELDI, 393 X NEWAVE,NORTH,WEST,ZMISS) 394 IF( IRET.NE.0 ) THEN 395 CALL INTLOG(JP_ERROR,'INTWAVE2: Interpolation failed.',JPQUIET) 396 INTWAVE2 = JPROUTINE + 13 397 GOTO 900 398 ENDIF 399C 400 ENDIF 401C 402C ------------------------------------------------------------------ 403C Section 5. Pack the interpolated field into GRIB. 404C ------------------------------------------------------------------ 405C 406 500 CONTINUE 407 408C If a 'bitmap' has been specified, build the bitmap 409C 410 LBITMP = LNOBITMP.AND. 411 X ((NOREPR.EQ.JPREGROT).OR.(NOREPR.EQ.JPREGULAR).OR. 412 X (NOREPR.EQ.JPGAUSSIAN)) 413 414 IF( LBITMP ) THEN 415 CALL INTLOG(JP_DEBUG,'INTWAVE2: MKBITMP is enabled',JPQUIET) 416 IERR = MKBITMP(NUM_E_W,NUM_N_S,NEWAVE,ZMISS) 417 IF( IERR.NE.0 ) THEN 418 CALL INTLOG(JP_ERROR,'INTFB: Problem applying bitmap',JPQUIET) 419 GOTO 900 420 ENDIF 421 ENDIF 422C If a 'frame' has been specified, build the frame 423C 424 LFRAME = LNOFRAME.AND. 425 X ((NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPGAUSSIAN).OR. 426 X (NOREPR.EQ.JPREGROT ).OR.(NOREPR.EQ.JPFGGROT ) ) 427 428 CALL INTLOG(JP_DEBUG,'INTWAVE2: NUM_E_W ',NUM_E_W) 429 CALL INTLOG(JP_DEBUG,'INTWAVE2: NUM_N_S ',NUM_N_S) 430 IF( LFRAME ) THEN 431 CALL INTLOG(JP_DEBUG,'INTWAVE2: MKFRAME is enabled',JPQUIET) 432 CALL INTLOG(JP_DEBUG,'INTWAVE2: NOFRAME ',NOFRAME) 433 CALL MKFRAME(NUM_E_W,NUM_N_S,NEWAVE,ZMISS,NOFRAME) 434 ENDIF 435C 436C GRIB sections 2 and 3 (bitmap) included 437C 438 ISEC1(5) = 192 439C 440 IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPREDLL) ) THEN 441C 442C Reduced gaussian or reduced lat/long with (possible) subarea 443C 444 IF( NOREPR.EQ.JPQUASI ) THEN 445 ISEC2(1) = JPGAUSSIAN 446 ISEC2(10)= NOGAUSS 447 ELSE 448 ISEC2(1) = JPREGULAR 449 ISEC2(10)= NINT(GRIDNS*1000.0) 450 ENDIF 451 ISEC2(2) = 255 452 ISEC2(3) = NUM_N_S 453 ISEC2(5) = NINT(WEST*1000.0) 454 ISEC2(6) = 0 455 ISEC2(8) = NINT(EAST*1000.0) 456 ISEC2(17)= 1 457 ISEC4(1) = 0 458C 459 LAT1 = 0 460 510 CONTINUE 461 LAT1 = LAT1 + 1 462 IF( NOREPR.EQ.JPQUASI ) THEN 463 RNLAT = ROGAUSS(LAT1) 464 RSLAT = ROGAUSS(LAT1+NUM_N_S-1) 465 ELSE 466 RNLAT = 90.0 - (LAT1-1)*GRIDNS 467 RSLAT = 90.0 - (LAT1+NUM_N_S-2)*GRIDNS 468 ENDIF 469 IF( RNLAT.GT.NORTH ) GOTO 510 470C 471 ISEC2(4) = NINT(RNLAT*1000.0) 472 ISEC2(7) = NINT(RSLAT*1000.0) 473C 474 DO LOOP = LAT1, (NUM_N_S+LAT1-1) 475 ISEC2(23+LOOP-LAT1) = NOLPTS(LOOP) 476 ISEC4(1) = ISEC4(1) + ISEC2(23+LOOP-LAT1) 477 ENDDO 478C 479 ELSE 480C 481C Regular gaussian or lat/long with (possible) subarea 482C 483 IF( NOREPR.EQ.JPGAUSSIAN ) THEN 484 ISEC2(1) = JPGAUSSIAN 485 ELSE 486 ISEC2(1) = JPREGULAR 487 ENDIF 488 ISEC2(2) = NUM_E_W 489 ISEC2(3) = NUM_N_S 490 ISEC2(4) = NINT(NORTH*1000.0) 491 ISEC2(5) = NINT(WEST*1000.0) 492 ISEC2(6) = 128 493 ISEC2(7) = NINT(SOUTH*1000.0) 494 ISEC2(8) = NINT(EAST*1000.0) 495 ISEC2(9) = NINT(GRIDWE*1000.0) 496 IF( NOREPR.EQ.JPGAUSSIAN ) THEN 497 ISEC2(10) = NOGAUSS 498 ELSE 499 ISEC2(10) = NINT(GRIDNS*1000.0) 500 ENDIF 501 ISEC2(17)= 0 502 ISEC4(1) = NUM_E_W * NUM_N_S 503 ENDIF 504C 505C 506 ISEC3(2) = NINT(ZMISS) 507 ZSEC3(2) = ZMISS 508C 509 ISEC4(1) = ISEC4(1) * NSPEC 510 ISEC4(2) = NOACC 511 IF( LCOEFFS ) ISEC4(6) = 16 512C `-----> wave fields have additional flags for 513C NC1 and NC2 coefficients (floats stored 514C in integer array isec4) 515 IRET = 1 516C 517C If grid-point output, setup for 2nd order packing if requested. 518C 519 IF( (NOREPR.NE.JPSPHERE) .AND. (NOREPR.NE.JPSPHROT) ) THEN 520 HFUNC = 'C' 521 IF( NOHFUNC.EQ.'K' ) THEN 522 HFUNC = 'K' 523 ISEC4(4) = 64 524 ISEC4(6) = 16 525 ISEC4(9) = 32 526 ISEC4(10) = 16 527 ISEC4(12) = 8 528 ISEC4(13) = 4 529 ISEC4(14) = 0 530 ISEC4(15) = -1 531 ELSE 532 ISEC4(4) = 0 533 ENDIF 534 ELSE 535 HFUNC = 'C' 536 IF( NOHFUNC.EQ.'C' ) THEN 537 ISEC2(6) = 2 538 ISEC4(4) = 64 539 ELSE IF( NOHFUNC.EQ.'S' ) THEN 540 ISEC2(6) = 1 541 ISEC4(4) = 0 542 ENDIF 543 544 ENDIF 545 546C 547 CALL GRIBEX (ISEC0,ISEC1,ISEC2,ZSEC2,ISEC3,ZSEC3,ISEC4, 548 X NEWAVE,ISEC4(1),OUTGRIB,OUTLEN,IWORD,HFUNC,IRET) 549C 550 IF ( IRET.NE.0 ) THEN 551 CALL INTLOG(JP_ERROR,'INTWAVE2: GRIBEX encoding failed.',IRET) 552 INTWAVE2 = JPROUTINE + 14 553 GOTO 900 554 ENDIF 555 OUTLEN = IWORD 556C 557C ------------------------------------------------------------------ 558C* Section 9. Closedown. 559C ------------------------------------------------------------------ 560C 561 900 CONTINUE 562C 563C Clear change flags for next product processing 564C 565 LCHANGE = .FALSE. 566 LSMCHNG = .FALSE. 567C 568 RETURN 569 END 570