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 HNTFAPH(INGRIB,INLEN) 12C 13C----> 14C**** HNTFAPH 15C 16C Purpose 17C ------- 18C 19C Prepare to interpolate to grid point field. 20C 21C 22C Interface 23C --------- 24C 25C IRET = HNTFAPH(INGRIB,INLEN) 26C 27C Input 28C ----- 29C 30C INGRIB - Input field (packed). 31C INLEN - Input field length (words). 32C 33C 34C Output 35C ------ 36C 37C Field unpacked values are in ZNFELDI, rotated if necessary. 38C 39C Returns: 0, if OK. Otherwise, an error occured in interpolation. 40C 41C 42C Method 43C ------ 44C 45C Unpack field. 46C 47C If the input is a spectral field and the output is a rotated 48C grid-point field, create a global reduced gaussian field and 49C then create the rotated grid-point field from it. 50C 51C 52C Externals 53C --------- 54C 55C GRIBEX - Decode/encode GRIB product. 56C GRSVCK - Turn off GRIB checking 57C INTLOG - Log error message. 58C INTLOGR - Log error message. 59C JDEBUG - Checks environment to switch on/off debug 60C FIXAREA - Fixup input/output field area definitions. 61C HSH2GG - Find suitable gaussian grid/spectral truncation 62C HIRLAM - Creates rotated lat/long field from reduced gaussian 63C HIRLSM - Creates rotated lat/long field from reduced gaussian 64C using land-sea mask 65C HRG2GG - Creates rotated gaussian field from reduced gaussian 66C HLL2LL - Creates rotated lat/long field from lat/long field 67C LSMFLD - Determines whether a field is to be interpolated using 68C a land-sea mask 69C 70C 71C Author 72C ------ 73C 74C J.D.Chambers ECMWF January 31, 2001 75C 76C----< 77C 78C -----------------------------------------------------------------| 79C* Section 0. Variables 80C -----------------------------------------------------------------| 81C 82 IMPLICIT NONE 83C 84C Function arguments 85C 86 INTEGER INGRIB(*),INLEN 87C 88#include "parim.h" 89#include "nifld.common" 90#include "nofld.common" 91#include "intf.h" 92#include "current.h" 93#include "grfixed.h" 94C 95C Parameters 96C 97 INTEGER JPROUTINE 98 PARAMETER (JPROUTINE = 40110 ) 99C 100C Local variables 101C 102 CHARACTER*1 HOLDTYP 103 CHARACTER*1 HTYPE 104 INTEGER IERR 105 INTEGER IRET 106 INTEGER ISHIZE 107 INTEGER ISIZE 108 INTEGER ITEMP,K,I,J,IM,JM 109 INTEGER IWORD 110 INTEGER KPTS(JPGTRUNC*2) 111 INTEGER LOOP 112 INTEGER NCOUNT 113 INTEGER NGAUSS 114 INTEGER NLAT 115 INTEGER NLON 116 INTEGER NTRUNC 117 INTEGER NUMPTS 118 INTEGER NUVFLAG 119 LOGICAL LLATOUT 120 LOGICAL LSP2RGG 121 LOGICAL LUSELSM 122 REAL AREA(4) 123 REAL EAST 124 REAL EW 125 REAL GRID(2) 126 REAL NS 127 REAL OLDGRID(2) 128 REAL PLATS(JPGTRUNC*2) 129 REAL POLE(2) 130 REAL RGGRID(1) 131 REAL SWORK(1) 132 REAL TEMP(1440,1440) 133 REAL WEST 134 REAL ZNFLDO(1) 135#ifndef _CRAYFTN 136#ifdef POINTER_64 137 INTEGER*8 IRGGRID 138 INTEGER*8 ISWORK 139 INTEGER*8 IZNFLDO 140#endif 141#endif 142 POINTER (IRGGRID, RGGRID) 143 POINTER (ISWORK, SWORK) 144 POINTER (IZNFLDO, ZNFLDO) 145C 146 SAVE IRGGRID 147 SAVE ISWORK 148 SAVE IZNFLDO 149C 150C Externals 151C 152 CHARACTER*1 GGHTYPE 153 INTEGER HSH2GG, HIRLAM, HIRLSM, HRG2GG, FIXAREA, PDDEFS 154 INTEGER HLL2LL 155 LOGICAL LSMFLD 156 EXTERNAL GGHTYPE 157 EXTERNAL HSH2GG, HIRLAM, HIRLSM, HRG2GG, FIXAREA, PDDEFS 158 EXTERNAL HLL2LL 159 EXTERNAL LSMFLD 160C 161C -----------------------------------------------------------------| 162C* Section 1. Initialise 163C -----------------------------------------------------------------| 164C 165 100 CONTINUE 166C 167 HNTFAPH = 0 168 LSP2RGG = .FALSE. 169C 170C -----------------------------------------------------------------| 171C* Section 2. Decode data from the GRIB code 172C -----------------------------------------------------------------| 173C 174 200 CONTINUE 175C 176C Decode data from GRIB code (no checking) 177C 178 IERR = 0 179 CALL GRSVCK(0) 180C 181 IERR = 1 182 IWORD = INLEN 183 ISEC3(2) = NINT(RMISSGV) 184 ZSEC3(2) = RMISSGV 185 186 CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4, 187 X ZNFELDI, JPEXPAND, INGRIB, INLEN, IWORD, 'D',IERR) 188 189 IF (ISEC2(11).EQ.64) THEN 190 CALL INTLOG(JP_DEBUG, 191 X 'HNTFAPH: Scanning flag west-east/south-north',ISEC2(11)) 192 ITEMP = NIAREA(1) 193 NIAREA(1) = NIAREA(3) 194 NIAREA(3) = ITEMP 195 196 IM = ISEC2(2) 197 JM = ISEC2(3) 198 K=0 199 DO J=JM,1,-1 200 DO I=1,IM 201 K=K+1 202 TEMP(I,J) = ZNFELDI(K) 203 END DO 204 END DO 205 K=0 206 DO J=1,JM 207 DO I=1,IM 208 K=K+1 209 ZNFELDI(K) = TEMP(I,J) 210 END DO 211 END DO 212 213 ENDIF 214 215C 216 IF( IERR.LT.0) THEN 217 IF( (IERR.EQ.-2).OR.(IERR.EQ.-4) ) THEN 218 CALL INTLOG(JP_DEBUG,'HNTFAPH: Use missing value',JPQUIET) 219 LIMISSV = .TRUE. 220 ELSE 221 CALL INTLOG(JP_ERROR,'HNTFAPH: GRIBEX decoding fail.',IERR) 222 HNTFAPH = JPROUTINE + 2 223 GOTO 900 224 ENDIF 225 ELSE IF( IERR.GT.0 ) THEN 226 CALL INTLOG(JP_ERROR,'HNTFAPH: GRIBEX decoding failed.',IERR) 227 HNTFAPH = JPROUTINE + 2 228 GOTO 900 229 ENDIF 230C 231 NCOUNT = ISEC4(1) 232C 233 LLATOUT = (NOREPR.EQ.JPREGROT).OR.(NOREPR.EQ.JPREGULAR) 234 235 IF( .NOT.LNOROTA ) GOTO 900 236C 237C -----------------------------------------------------------------| 238C* Section 3. Handle rotation, if necessary. 239C -----------------------------------------------------------------| 240C 241 300 CONTINUE 242C 243 CALL INTLOG(JP_DEBUG,'HNTFAPH: Rotate field.',JPQUIET) 244 CALL INTLOG(JP_DEBUG,'HNTFAPH: South pole lat ',NOROTA(1)) 245 CALL INTLOG(JP_DEBUG,'HNTFAPH: South pole long ',NOROTA(2)) 246C 247C Fill area limits (handles case when default 0/0/0/0 given) 248C 249 IRET = FIXAREA() 250 IF( IRET.NE.0 ) THEN 251 CALL INTLOG(JP_ERROR,'HNTFAPH: area fixup failed',JPQUIET) 252 HNTFAPH = JPROUTINE + 3 253 GOTO 900 254 ENDIF 255C 256 AREA(1) = REAL(NOAREA(1))/PPMULT 257 AREA(2) = REAL(NOAREA(2))/PPMULT 258 AREA(3) = REAL(NOAREA(3))/PPMULT 259 AREA(4) = REAL(NOAREA(4))/PPMULT 260C 261 GRID(1) = REAL(NOGRID(1))/PPMULT 262 GRID(2) = REAL(NOGRID(2))/PPMULT 263C 264 POLE(1) = REAL(NOROTA(1))/PPMULT 265 POLE(2) = REAL(NOROTA(2))/PPMULT 266C 267C -----------------------------------------------------------------| 268C* Section 4. Spectral to rotated grid-point 269C -----------------------------------------------------------------| 270C 271 400 CONTINUE 272C 273 IF( (NIREPR.EQ.JPSPHERE).OR.(NIREPR.EQ.JPSPHROT) ) THEN 274C 275C Convert spectral to suitable global reduced gaussian grid 276C 277 CALL INTLOG(JP_DEBUG, 278 X 'HNTFAPH: Spectral to suitable reduced gaussian',JPQUIET) 279C 280 NIRESO = ISEC2(2) 281 NTRUNC = ISEC2(2) 282 IF( LNORESO ) THEN 283 NTRUNC = NORESO 284 NGAUSS = 0 285 HTYPE = '' 286 NS = 0. 287 EW = 0. 288 IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE) 289 IF( IRET.NE.0 ) THEN 290 CALL INTLOG(JP_ERROR, 291 X 'HNTFAPH: problem getting data for reduced grid',NTRUNC) 292 HNTFAPH = JPROUTINE + 4 293 GOTO 900 294 ENDIF 295 GOTO 401 296 ENDIF 297 298 IF( LARESOL ) THEN 299 NS = 0. 300 EW = 0. 301 IF( LLATOUT ) THEN 302 NS = GRID(1) 303 EW = GRID(2) 304 NTRUNC = 0 305 NGAUSS = 0 306 HTYPE = '' 307 ELSE 308 HTYPE = 'R' 309 ENDIF 310 IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE) 311 IF( IRET.NE.0 ) THEN 312 CALL INTLOG(JP_ERROR, 313 X 'HNTFAPH: problem getting data for reduced grid',NTRUNC) 314 HNTFAPH = JPROUTINE + 4 315 GOTO 900 316 ENDIF 317 ENDIF 318 319C Truncate if a smaller resolution has been requested 320C 321 401 CONTINUE 322 IF( NTRUNC.LT.NIRESO ) THEN 323 CALL INTLOG(JP_DEBUG,'HNTFAPH: Truncation changed from:',NIRESO) 324 CALL INTLOG(JP_DEBUG,'HNTFAPH: to: ',NTRUNC) 325 CALL INTLOG(JP_DEBUG,'HNTFAPH: Gaussian number is:',NGAUSS) 326C 327 ISHIZE = (NTRUNC+1)*(NTRUNC+4) 328 CALL JMEMHAN( 3, IZNFLDO, ISHIZE, 1, IERR) 329 IF( IERR.NE.0 ) THEN 330 CALL INTLOG(JP_FATAL, 331 X 'HNTFAPH: Get scratch space failed',JPQUIET) 332 HNTFAPH = JPROUTINE + 4 333 GOTO 900 334 ENDIF 335 336C Generate spherical harmonics with output truncation 337 CALL SH2SH( ZNFELDI, NIRESO, ZNFLDO, NTRUNC ) 338 339C Move new spherical harmonics to 'input' array 340 ZNFELDI(1:ISHIZE) = ZNFLDO(1:ISHIZE) 341 342 ELSE 343 CALL INTLOG(JP_DEBUG, 344 X 'HNTFAPH: Spectral to suitable reduced gaussian',JPQUIET) 345C 346 NTRUNC = ISEC2(2) 347 NGAUSS = 0 348 HTYPE = '' 349 NS = 0. 350 EW = 0. 351 IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE) 352 IF( IRET.NE.0 ) THEN 353 CALL INTLOG(JP_ERROR, 354 X 'HNTFAPH: problem getting data for reduced grid',NTRUNC) 355 HNTFAPH = JPROUTINE + 4 356 GOTO 900 357 ENDIF 358 359 ENDIF 360 361C 362C Dynamically allocate memory for global reduced gaussian grid 363C 364 CALL JMEMHAN( 18, IRGGRID, ISIZE, 1, IRET) 365 IF( IRET.NE.0 ) THEN 366 CALL INTLOG(JP_ERROR, 367 X 'HNTFAPH: memory alloc for reduced grid fail',JPQUIET) 368 HNTFAPH = JPROUTINE + 4 369 GOTO 900 370 ENDIF 371 372C Set flag to show field is not a wind component 373 NUVFLAG = 0 374 375C Create the reduced gaussian grid 376 HOLDTYP = HOGAUST 377 WEST = 0.0 378 EAST = 360.0 - (360.0/(NGAUSS*4)) 379c EMOS-199: adjusted for reduced_gg/octahedral 380 IF (HTYPE.EQ.'O') EAST = 360.0 - (360.0/FLOAT(KPTS(NGAUSS))) 381 IF (HTYPE.NE.'R' .AND. HTYPE.NE.'O' .AND. 382 X HTYPE.NE.'F' .AND. HTYPE.NE.'U') THEN 383 HTYPE = 'R' 384 ENDIF 385 CALL JAGGGP(ZNFELDI,NTRUNC,PLATS(1),PLATS(NGAUSS*2),WEST, 386 X EAST,NGAUSS,HTYPE,KPTS,RGGRID,NUVFLAG,IRET) 387 IF( IRET.NE.0 ) THEN 388 CALL INTLOG(JP_ERROR, 389 X 'HNTFAPH: spectral to reduced gaussian failed',JPQUIET) 390 HNTFAPH = JPROUTINE + 4 391 GOTO 900 392 ENDIF 393 HOGAUST = HOLDTYP 394 395 NCOUNT = 0 396 DO LOOP = 1, (NGAUSS*2) 397 NCOUNT= NCOUNT + KPTS(LOOP) 398 ENDDO 399 400 LSP2RGG = .TRUE. 401 402 IF( NOREPR.EQ.JPFGGROT ) THEN 403 CALL INTLOG(JP_DEBUG, 404 X 'HNTFAPH: Convert gaussian to rotated gaussian',JPQUIET) 405 GOTO 600 406 ENDIF 407 408 ENDIF 409C 410C -----------------------------------------------------------------| 411C* Section 5. Gaussian to rotated lat/long 412C -----------------------------------------------------------------| 413C 414 500 CONTINUE 415C 416 IF( (LLATOUT) .AND. ( 417 X (NIREPR.EQ.JPQUASI) .OR. 418 X (NIREPR.EQ.JPGAUSSIAN) .OR. 419 X (LSP2RGG) )) THEN 420C 421 CALL INTLOG(JP_DEBUG,'HNTFAPH: Gauss to lat/lon',JPQUIET) 422C 423C Dynamically allocate memory for rotated lat/long grid 424C 425 NLON = 1 + NINT((AREA(JPEAST) - AREA(JPWEST)) / 426 X GRID(JPWESTEP)) !SC 427 NLAT = 1 + NINT((AREA(JPNORTH) - AREA(JPSOUTH)) / 428 X GRID(JPNSSTEP)) !SC 429C 430 NOWE = NLON 431 NONS = NLAT 432C 433 NUMPTS = NLON * NLAT 434 ISIZE = NUMPTS 435 CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET) 436 IF( IRET.NE.0 ) THEN 437 CALL INTLOG(JP_ERROR, 438 X 'HNTFAPH: memory alloc for lat/long grid fail',JPQUIET) 439 HNTFAPH = JPROUTINE + 5 440 GOTO 900 441 ENDIF 442C 443 LUSELSM = LSMFLD() 444C 445C If original field was spectral, ... 446C 447 IF( LSP2RGG ) THEN 448 IF( LUSELSM ) THEN 449 IRET = HIRLSM(LO12PT,RGGRID,NCOUNT,NGAUSS,HTYPE,AREA,POLE, 450 X GRID,SWORK,ISIZE,NLON,NLAT) 451 ELSE 452 IRET = HIRLAM(LO12PT,RGGRID,NCOUNT,NGAUSS,HTYPE,AREA,POLE, 453 X GRID,SWORK,ISIZE,NLON,NLAT) 454 ENDIF 455C 456 ELSE 457C 458C If original field was gaussian, ... 459C 460 IRET = PDDEFS() 461 NGAUSS = ISEC2(10) 462 HTYPE = GGHTYPE(NIREPR,NGAUSS,MILLEN) 463cs IF( LUSELSM ) THEN 464 IF( LSM ) THEN 465 IRET = HIRLSM(LO12PT,ZNFELDI,NCOUNT,NGAUSS,HTYPE,AREA,POLE, 466 X GRID,SWORK,ISIZE,NLON,NLAT) 467 ELSE 468 IRET = HIRLAM(LO12PT,ZNFELDI,NCOUNT,NGAUSS,HTYPE,AREA,POLE, 469 X GRID,SWORK,ISIZE,NLON,NLAT) 470 ENDIF 471C 472 ENDIF 473C 474 IF( IRET.NE.0 ) THEN 475 CALL INTLOG(JP_ERROR, 476 X 'HNTFAPH: HIRLAM rotation failed',JPQUIET) 477 HNTFAPH = JPROUTINE + 7 478 GOTO 900 479 ENDIF 480C 481 ISEC2(1) = JPREGROT 482 ISEC4(1) = NOWE * NONS 483C 484 ENDIF 485C 486C -----------------------------------------------------------------| 487C* Section 6. Gaussian to rotated gaussian 488C -----------------------------------------------------------------| 489C 490 600 CONTINUE 491C 492 IF( (LSP2RGG.AND.(NOREPR.EQ.JPFGGROT)) .OR. 493 X (LSP2RGG.AND.(NOREPR.EQ.JPQGGROT)) .OR. 494 X (((NIREPR.EQ.JPQUASI).OR.(NIREPR.EQ.JPGAUSSIAN)) .AND. 495 X (NOREPR.EQ.JPFGGROT).OR.(NOREPR.EQ.JPQGGROT)) ) THEN 496 CALL INTLOG(JP_DEBUG, 497 X 'HNTFAPH: Gaussian to reduced gaussian',JPQUIET) 498C 499C Dynamically allocate memory for rotated lat/long grid 500C 501C ISIZE = NOGAUSS * NOGAUSS * 8 502 ISIZE = (2*NOGAUSS) * (4*NOGAUSS + 20) ! account for RGG/octahedral 503 CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET) 504 IF( IRET.NE.0 ) THEN 505 CALL INTLOG(JP_ERROR, 506 X 'HNTFAPH: memory alloc for gaussian grid fail',JPQUIET) 507 HNTFAPH = JPROUTINE + 6 508 GOTO 900 509 ENDIF 510C 511C If original field was spectral, ... 512C 513 IF( LSP2RGG ) THEN 514 IRET = HRG2GG(LO12PT,RGGRID,NGAUSS,AREA,POLE, 515 X NOGAUSS,HOGAUST,SWORK,ISIZE,NUMPTS) 516C 517 ELSE 518C 519C If original field was gaussian, ... 520C 521 NGAUSS = ISEC2(10) 522 IRET = HRG2GG(LO12PT,ZNFELDI,NGAUSS,AREA,POLE, 523 X NOGAUSS,HOGAUST,SWORK,ISIZE,NUMPTS) 524 ENDIF 525 IF( IRET.NE.0 ) THEN 526 CALL INTLOG(JP_ERROR, 527 X 'HNTFAPH: HRG2GG rotation failed',JPQUIET) 528 HNTFAPH = JPROUTINE + 6 529 GOTO 900 530 ENDIF 531C 532 IF( (NOREPR.EQ.JPQUASI).OR. 533 X (NOREPR.EQ.JPQGGROT).OR. 534 X (NOREPR.EQ.JPFGGROT).OR. 535 X (NOREPR.EQ.JPGAUSSIAN) ) THEN 536 ISEC2(1) = JPFGGROT 537 ELSE 538 ISEC2(1) = NOREPR 539 ENDIF 540C 541 ISEC4(1) = NUMPTS 542C 543 ENDIF 544C 545C -----------------------------------------------------------------| 546C* Section 7. Lat/long to rotated lat/long 547C -----------------------------------------------------------------| 548C 549 700 CONTINUE 550C 551 IF( (NIREPR.EQ.JPREGULAR) ) THEN 552C 553C Dynamically allocate memory for rotated lat/long grid 554C 555 NLON = 1 + NINT((AREA(JPEAST) - AREA(JPWEST)) / 556 X GRID(JPWESTEP)) 557 NLAT = 1 + NINT((AREA(JPNORTH) - AREA(JPSOUTH)) / 558 X GRID(JPNSSTEP)) 559C 560 NOWE = NLON 561 NONS = NLAT 562C 563 NUMPTS = NLON * NLAT 564 ISIZE = NUMPTS 565 CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET) 566 IF( IRET.NE.0 ) THEN 567 CALL INTLOG(JP_ERROR, 568 X 'HNTFAPH: memory alloc for lat/long grid fail',JPQUIET) 569 HNTFAPH = JPROUTINE + 7 570 GOTO 900 571 ENDIF 572C 573 OLDGRID(1) = REAL(NIGRID(1))/100000.0 574 OLDGRID(2) = REAL(NIGRID(2))/100000.0 575 IRET = HLL2LL(LO12PT,ZNFELDI,OLDGRID,AREA,POLE,GRID,SWORK,ISIZE, 576 X NLON,NLAT) 577C 578 IF( IRET.NE.0 ) THEN 579 CALL INTLOG(JP_ERROR, 580 X 'HNTFAPH: HLL2LL rotation failed',JPQUIET) 581 HNTFAPH = JPROUTINE + 7 582 GOTO 900 583 ENDIF 584C 585 ISEC2(1) = JPREGROT 586 ISEC4(1) = NOWE * NONS 587C 588 ENDIF 589 590C move rotated field back into field original array 591 ZNFELDI(1:NUMPTS) = SWORK(1:NUMPTS) 592 593C -----------------------------------------------------------------| 594C* Section 9. Closedown. 595C -----------------------------------------------------------------| 596C 597 900 CONTINUE 598C 599 RETURN 600 END 601