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 HRG2GG(L12PNT,OLDFLD,KGNOLD,AREA,POLE, 12 X KGNNEW,HTYPE,NEWFLD,KSIZE,NUMPTS) 13C 14C----> 15C**** HRG2GG 16C 17C Purpose 18C ------- 19C 20C This routine creates a rotated gaussian field from a reduced 21C gaussian field using 12-point horizontal interpolation. 22C 23C 24C Interface 25C --------- 26C 27C IRET = HRG2GG(L12PNT,OLDFLD,KGNOLD,AREA,POLE, 28C X KGNNEW,HTYPE,NEWFLD,KSIZE,NUMPTS) 29C 30C 31C Input parameters 32C ---------------- 33C 34C L12PNT - Chooses between 12-point and 4-point interpolation 35C OLDFLD - Array of values from the old (reduced) gaussian field 36C KGNOLD - Gaussian number for the old gaussian field 37C AREA - Limits of area (N/W/S/E) 38C POLE - Pole of rotation (lat/long) 39C KGNNEW - Gaussian grid number (for the new field) 40C HTYPE - Gaussian grid type (for the new field) 41C = 'R' for reduced ("quasi-regular"), 42C = 'O' for reduced/octahedral, 43C = 'F' for full 44C KSIZE - The size of the array to fill with the regular 45C lat/long field 46C 47C 48C Output parameters 49C ----------------- 50C 51C NEWFLD - The array of values for the gaussian field 52C NUMPTS - Number of points in the new gaussian field. 53C 54C Returns 0 if function successful, non-zero otherwise. 55C 56C 57C Common block usage 58C ------------------ 59C 60C None 61C 62C 63C Method 64C ------ 65C 66C Numbering of the points (I is the interpolation point): 67C 68C 13 5 6 14 69C 70C 7 1 2 8 71C (I) 72C 9 3 4 10 73C 74C 15 11 12 16 75C 76C The 12-point interpolation is not possible if either of the top 77C two rows is above the original field northern latitude. The 78C nearest neighbour is used if both rows are above, and a 4-pt 79C bilinear interpolation is used if the top row is above. 80C Similarily, if either of the bottom two rows is below the original 81C field southern latitude. 82C 83C 84C Externals 85C --------- 86C 87C INTLOG - Logs messages 88C JMALLOC - Dynamically allocate memory 89C JFREE - Free dynamically allocated memory 90C JGETGG - Reads the definition of a gaussian grid 91C HGENGG - Calculates original lat/long (before rotation) for 92C a rotated gaussian grid 93C HNEI12 - Finds neighbours for points for interpolation 94C HWTS12 - Calculates weightings for points for interpolation 95C FORCED_NEAREST_NEIGHBOUR - check forced interpolation method 96C 97C 98C Reference 99C --------- 100C 101C None. 102C 103C 104C Comments 105C -------- 106C 107C None. 108C 109C 110C Author 111C ------ 112C 113C J.D.Chambers ECMWF February 2001 114C 115C 116C Modifications 117C ------------- 118C 119C None. 120C 121C----< 122C -----------------------------------------------------------------| 123C* Section 0. Definition of variables. 124C -----------------------------------------------------------------| 125C 126 IMPLICIT NONE 127C 128#include "parim.h" 129#include "nifld.common" 130#include "nofld.common" 131C 132C Parameters 133C 134 INTEGER JNORTH, JSOUTH, JWEST, JEAST, JW_E, JN_S, JLAT, JLON 135 INTEGER JP12PT, JP4PT, JPNEARN 136 PARAMETER (JP12PT = 0) 137 PARAMETER (JP4PT = 1) 138 PARAMETER (JPNEARN = 2) 139 PARAMETER (JNORTH = 1 ) 140 PARAMETER (JWEST = 2 ) 141 PARAMETER (JSOUTH = 3 ) 142 PARAMETER (JEAST = 4 ) 143 PARAMETER (JW_E = 1 ) 144 PARAMETER (JN_S = 2 ) 145 PARAMETER (JLAT = 1 ) 146 PARAMETER (JLON = 2 ) 147C 148C Function arguments 149C 150 LOGICAL L12PNT 151 INTEGER KGNOLD, KGNNEW, KSIZE, NUMPTS 152 REAL AREA(4), POLE(2), OLDFLD(*), NEWFLD(KSIZE) 153 CHARACTER*1 HTYPE 154C 155C Local variables 156C 157 INTEGER NEXT, LOOP, IRET, NLEN, NOPREV, NNPREV, NBYTES, NUMBER 158 INTEGER NEAREST 159 CHARACTER*1 NNTYPE 160C 161 LOGICAL LNEW, LFIRST, LVEGGY 162 INTEGER KSCHEME(1),NEIGH(12,1), KLA(1) 163 REAL PWTS(12,1) 164 POINTER (IPKSCHE, KSCHEME) 165 POINTER (IPNEIGH, NEIGH) 166 POINTER (IPKLA, KLA) 167 POINTER (IPPWTS, PWTS) 168C 169 REAL PDLO0(1),PDLO1(1),PDLO2(1),PDLO3(1),PDLAT(1) 170 POINTER (IPPDLO0, PDLO0) 171 POINTER (IPPDLO1, PDLO1) 172 POINTER (IPPDLO2, PDLO2) 173 POINTER (IPPDLO3, PDLO3) 174 POINTER (IPPDLAT, PDLAT) 175C 176 INTEGER IGG, IGGOLD, INN, INNOLD 177 INTEGER KPTS(1), KNPTS(1) 178 REAL GLATS(1), GNLATS(1) 179 INTEGER IOFFS(1) 180 POINTER (IPKPTS, KPTS) 181 POINTER (IPKNPTS, KNPTS) 182 POINTER (IPIOFFS, IOFFS) 183 POINTER (IPGLATS, GLATS) 184 POINTER (IPGNLATS, GNLATS) 185C 186 INTEGER ILL, ILLOLD 187 REAL RLAT(1),RLON(1) 188 POINTER (IPRLAT, RLAT) 189 POINTER (IPRLON, RLON) 190C 191 REAL OLD(1) 192 POINTER (IOLD, OLD) 193C 194 DATA NOPREV/-1/, NNPREV/-1/, NNTYPE/' '/ 195 DATA LNEW/.FALSE./, LFIRST/.TRUE./ 196 DATA IGGOLD/-1/, INNOLD/-1/, ILLOLD/-1/, IOLD/-1/ 197C 198 SAVE LNEW, LFIRST 199 SAVE IPKSCHE, IPNEIGH, IPKLA, IPPWTS 200 SAVE IPPDLO0, IPPDLO1, IPPDLO2, IPPDLO3, IPPDLAT 201 SAVE NOPREV,NNPREV,NNTYPE 202 SAVE IGGOLD,INNOLD,IPKPTS,IPKNPTS,IPIOFFS,IPGLATS,IPGNLATS 203 SAVE ILLOLD, IPRLAT, IPRLON, IOLD 204C 205C Externals 206C 207 LOGICAL FORCED_NEAREST_NEIGHBOUR 208 INTEGER HNEI12, HGENGG 209#ifdef POINTER_64 210 INTEGER*8 JMALLOC 211#else 212 INTEGER JMALLOC 213#endif 214C 215C -----------------------------------------------------------------| 216C Section 1. Initialise. 217C -----------------------------------------------------------------| 218C 219 100 CONTINUE 220C 221 HRG2GG = 0 222C 223 CALL JDEBUG() 224 225 IF( L12PNT ) THEN 226 CALL INTLOG(JP_DEBUG,'HRG2GG: 12-pt interpolation',JPQUIET) 227 ELSE 228 CALL INTLOG(JP_DEBUG,'HRG2GG: 4-pt interpolation',JPQUIET) 229 ENDIF 230C 231 CALL CHKPREC() 232 IF( LPREC )THEN 233 CALL INTLOG(JP_DEBUG, 234 X 'HRG2GG: precipitation threshold applied',JPQUIET) 235 ELSE 236 CALL INTLOG(JP_DEBUG, 237 X 'HRG2GG: precipitation threshold not applied',JPQUIET) 238 ENDIF 239 240C Use nearest neighbour if required 241 LVEGGY = FORCED_NEAREST_NEIGHBOUR(LMETHOD,NITABLE,NIPARAM) 242 IF( LVEGGY ) CALL INTLOG(JP_DEBUG, 243 X 'HRG2GG: nearest neighbour processing',JPQUIET) 244 245 246C 247C Dynamically allocate memory for old gaussian grid information. 248C 249 IGG = KGNOLD*2 250C 251 IF( IGG.GT.IGGOLD ) THEN 252C 253 IF( IGGOLD.GT.0 ) CALL JFREE(IPKPTS) 254C 255 NBYTES = (IGG*JPRLEN) + (2*IGG+1)*JPBYTES 256C 257 IPKPTS = JMALLOC(NBYTES) 258#ifdef hpR64 259 IPKPTS = IPKPTS/(1024*1024*1024*4) 260#endif 261 IF( IPKPTS.EQ.0 ) THEN 262 CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET) 263 HRG2GG = 1 264 GOTO 900 265 ENDIF 266C 267 IPGLATS = IPKPTS + (IGG*JPBYTES) 268 IPIOFFS = IPGLATS + (IGG*JPRLEN) 269C 270 IGGOLD = IGG 271 NOPREV = -1 272C 273 ENDIF 274C 275C Build up offsets to start of each latitude in the original field. 276C 277 IF( KGNOLD.NE.NOPREV ) THEN 278 CALL JGETGG(KGNOLD,HTYPE,GLATS,KPTS,IRET) 279 IF( IRET.NE.0 ) THEN 280 CALL INTLOG(JP_ERROR, 281 X 'HRG2GG: JGETGG failed to get gaussian data',JPQUIET) 282 HRG2GG = 1 283 GOTO 900 284 ENDIF 285C 286 IOFFS(1) = 1 287 DO LOOP = 2, (KGNOLD*2+1) 288 IOFFS(LOOP) = IOFFS(LOOP-1) + KPTS(LOOP-1) 289 ENDDO 290C 291C Allocate memory to hold the input field 292C (in case OLDFLD and NEWFLD are the same arrays) 293C 294 IF( IOLD.GT.0 ) CALL JFREE(IOLD) 295C 296 NUMBER = (IOFFS(KGNOLD*2+1) - 1) 297 NBYTES = NUMBER * JPRLEN 298C 299 IOLD = JMALLOC(NBYTES) 300#ifdef hpR64 301 IOLD = IOLD/(1024*1024*1024*4) 302#endif 303 IF( IOLD.EQ.0 ) THEN 304 CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET) 305 HRG2GG = 1 306 GOTO 900 307 ENDIF 308C 309 NOPREV = KGNOLD 310 ENDIF 311C 312C Preserve the input field 313C 314 NUMBER = (IOFFS(KGNOLD*2+1) - 1) 315 DO LOOP = 1, NUMBER 316 OLD(LOOP) = OLDFLD(LOOP) 317 ENDDO 318C 319C -----------------------------------------------------------------| 320C Section 2. Generate the lat/long points for the output grid 321C -----------------------------------------------------------------| 322C 323 200 CONTINUE 324C 325C 326C Dynamically allocate memory for gaussian grid information. 327C 328 INN = KGNNEW*2 329C 330 IF( INN.GT.INNOLD ) THEN 331C 332 IF( INNOLD.GT.0 ) CALL JFREE(IPKNPTS) 333C 334 NBYTES = (INN*JPRLEN) + (INN*JPBYTES) 335C 336 IPKNPTS = JMALLOC(NBYTES) 337#ifdef hpR64 338 IPKNPTS = IPKNPTS/(1024*1024*1024*4) 339#endif 340 IF( IPKNPTS.EQ.0 ) THEN 341 CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET) 342 HRG2GG = 2 343 GOTO 900 344 ENDIF 345C 346 IPGNLATS = IPKNPTS + (INN*JPBYTES) 347C 348 INNOLD = INN 349 NNPREV = -1 350 ENDIF 351C 352 IF( (KGNNEW.NE.NNPREV).OR.(HTYPE.NE.NNTYPE) ) THEN 353 IF( HTYPE.NE.'F' .AND. HTYPE.NE.'f' .AND. 354 X HTYPE.NE.'O' .AND. HTYPE.NE.'o' ) THEN 355 HTYPE = 'R' 356 ENDIF 357 CALL JGETGG(KGNNEW,HTYPE,GNLATS,KNPTS,IRET) 358 IF( IRET.NE.0 ) THEN 359 CALL INTLOG(JP_ERROR, 360 X 'HRG2GG: JGETGG failed to get gaussian data',JPQUIET) 361 HRG2GG = 2 362 GOTO 900 363 ENDIF 364C 365 IF( HTYPE.EQ.'F' .OR. HTYPE.EQ.'f' ) THEN 366 NLEN = KGNNEW * KGNNEW * 8 367 ELSE 368 NLEN = 0 369 DO LOOP = 1, KGNNEW*2 370 NLEN = NLEN + KNPTS(LOOP) 371 ENDDO 372 ENDIF 373 NNPREV = KGNNEW 374 NNTYPE = HTYPE 375 ENDIF 376C 377C Check that given array is big enough for the new field. 378C 379 IF( NLEN.GT.KSIZE ) THEN 380 CALL INTLOG(JP_ERROR,'HRG2GG: Given array size = ',KSIZE) 381 CALL INTLOG(JP_ERROR,'HRG2GG: Required size = ',NLEN) 382 HRG2GG = 2 383 GOTO 900 384 ENDIF 385C 386C Dynamically allocate memory for lat/long arrays. 387C 388 ILL = NLEN 389 IF( ILL.GT.ILLOLD ) THEN 390C 391 LNEW = .TRUE. 392C 393 IF( ILLOLD.GT.0 ) CALL JFREE(IPRLON) 394C 395 NBYTES = 2*ILL*JPRLEN 396C 397 IPRLON = JMALLOC(NBYTES) 398#ifdef hpR64 399 IPRLON = IPRLON/(1024*1024*1024*4) 400#endif 401 IF( IPRLON.EQ.0 ) THEN 402 CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET) 403 HRG2GG = 2 404 GOTO 900 405 ENDIF 406C 407 IPRLAT = IPRLON + (ILL*JPRLEN) 408C 409 ILLOLD = ILL 410C 411 ENDIF 412C 413 IRET = HGENGG(AREA,POLE,KGNNEW,HTYPE,KNPTS,GNLATS,ILL, 414 X RLAT,RLON,NUMPTS) 415 IF( IRET.NE.0 ) THEN 416 CALL INTLOG(JP_ERROR, 417 X 'HRG2GG: HGENGG failed to get lat/lon grid data',JPQUIET) 418 HRG2GG = 2 419 GOTO 900 420 ENDIF 421C 422C -----------------------------------------------------------------| 423C Section 3. Find neighbours for each point for interpolation. 424C -----------------------------------------------------------------| 425C 426 300 CONTINUE 427C 428C Dynamically allocate memory for interpolation arrays. 429C 430 IF( LNEW ) THEN 431C 432 IF( .NOT.LFIRST ) CALL JFREE(IPPDLO0) 433C 434 NBYTES = (17*JPRLEN + 14*JPBYTES) * ILL 435C 436 IPPDLO0 = JMALLOC(NBYTES) 437#ifdef hpR64 438 IPPDLO0 = IPPDLO0/(1024*1024*1024*4) 439#endif 440 IF( IPPDLO0.EQ.0 ) THEN 441 CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET) 442 HRG2GG = 3 443 GOTO 900 444 ENDIF 445C 446 IPPDLO1 = IPPDLO0 + (ILL*JPRLEN) 447 IPPDLO2 = IPPDLO1 + (ILL*JPRLEN) 448 IPPDLO3 = IPPDLO2 + (ILL*JPRLEN) 449 IPPDLAT = IPPDLO3 + (ILL*JPRLEN) 450 IPPWTS = IPPDLAT + (ILL*JPRLEN) 451 IPKSCHE = IPPWTS + (12*ILL*JPRLEN) 452 IPKLA = IPKSCHE + (ILL*JPBYTES) 453 IPNEIGH = IPKLA + (ILL*JPBYTES) 454C 455 LFIRST = .FALSE. 456 LNEW = .FALSE. 457C 458 ENDIF 459C 460C Find neighbours. 461C 462 IRET = HNEI12(L12PNT,NLEN,RLAT,RLON,KGNOLD,KPTS,GLATS, 463 X KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH) 464 IF( IRET.NE.0 ) THEN 465 CALL INTLOG(JP_ERROR, 466 X 'HRG2GG: HNEI12 failed to find neighbours',JPQUIET) 467 HRG2GG = 3 468 GOTO 900 469 ENDIF 470C 471C -----------------------------------------------------------------| 472C Section 4. Perform the 12-point horizontal interpolation. 473C -----------------------------------------------------------------| 474C 475 400 CONTINUE 476C 477C Setup the 12-point horizontal interpolation weights 478C 479 CALL HWTS12 480 X (NLEN,KSCHEME,KLA,PDLAT,GLATS,PDLO0,PDLO1,PDLO2,PDLO3,NEIGH, 481 X PWTS) 482C 483C Calculate the interpolated grid point values 484C 485 DO LOOP = 1, NLEN 486 IF( LVEGGY) THEN 487 NEAREST = 1 488 IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2 489 IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3 490 IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4 491 IF( PWTS( 5,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 5 492 IF( PWTS( 6,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 6 493 IF( PWTS( 7,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 7 494 IF( PWTS( 8,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 8 495 IF( PWTS( 9,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 9 496 IF( PWTS(10,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =10 497 IF( PWTS(11,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =11 498 IF( PWTS(12,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =12 499 NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP)) 500 ELSE 501 IF( KSCHEME(LOOP).EQ.JP12PT ) THEN 502 NEWFLD(LOOP) = 503 X OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) + 504 X OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) + 505 X OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) + 506 X OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP) + 507 X OLD(NEIGH( 5,LOOP)) * PWTS( 5,LOOP) + 508 X OLD(NEIGH( 6,LOOP)) * PWTS( 6,LOOP) + 509 X OLD(NEIGH( 7,LOOP)) * PWTS( 7,LOOP) + 510 X OLD(NEIGH( 8,LOOP)) * PWTS( 8,LOOP) + 511 X OLD(NEIGH( 9,LOOP)) * PWTS( 9,LOOP) + 512 X OLD(NEIGH(10,LOOP)) * PWTS(10,LOOP) + 513 X OLD(NEIGH(11,LOOP)) * PWTS(11,LOOP) + 514 X OLD(NEIGH(12,LOOP)) * PWTS(12,LOOP) 515C 516 ELSE IF( KSCHEME(LOOP).EQ.JP4PT ) THEN 517 NEWFLD(LOOP) = 518 X OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) + 519 X OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) + 520 X OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) + 521 X OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP) 522C 523 ELSE 524 DO NEXT = 1, 4 525 IF( NEIGH(NEXT,LOOP).NE.0 ) 526 X NEWFLD(LOOP) = OLD(NEIGH(NEXT,LOOP)) 527 ENDDO 528C 529 ENDIF 530 ENDIF 531 ENDDO 532C 533C -----------------------------------------------------------------| 534C Section 9. Return. 535C -----------------------------------------------------------------| 536C 537 900 CONTINUE 538C 539 RETURN 540 END 541 542