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 HIRLAM(L12PNT,OLDFLD,KOUNT,KGAUSS, 12 X HTYPE,AREA,POLE,GRID,NEWFLD,KSIZE,NLON,NLAT) 13C 14C----> 15C**** HIRLAM 16C 17C Purpose 18C ------- 19C 20C This routine creates a rotated regular lat/long field from a 21C gaussian field using 12-point horizontal interpolation. 22C 23C 24C Interface 25C --------- 26C 27C IRET = HIRLAM(L12PNT,OLDFLD,KOUNT,KGAUSS,HTYPE,AREA,POLE,GRID, 28C X NEWFLD,KSIZE,NLON,NLAT) 29C 30C 31C Input parameters 32C ---------------- 33C 34C L12PNT - Chooses between 12-point and 4-point interpolation 35C = .TRUE. for 12-point horizontal 36C = .FALSE. for 4-point 37C OLDFLD - Array of values from the (old) gaussian field 38C KOUNT - Number of values in OLDFLD 39C KGAUSS - Gaussian grid number 40C HTYPE - Gaussian grid type 41C = 'R' for reduced ("quasi-regular"), 42C = 'O' for reduced/octahedral, 43C = 'F' for full, 44C = 'U' for a user-defined gaussian grid 45C AREA - Limits of output area (N/W/S/E) 46C POLE - Pole of rotation (lat/long) 47C GRID - Output lat/long grid increments (we/ns) 48C KSIZE - The size of the output array to fill with the regular 49C lat/long field 50C 51C 52C Output parameters 53C ----------------- 54C 55C NEWFLD - The array of values for the regular lat/long field 56C NLON - Number of longitudes in the regular lat/long field 57C NLAT - Number of latitudes in the regular lat/long field 58C 59C Returns 0 if function successful, non-zero otherwise. 60C 61C Common block usage 62C ------------------ 63C 64C nifld.common 65C nofld.common 66C 67C 68C Method 69C ------ 70C 71C Numbering of the points (I is the interpolation point): 72C 73C 13 5 6 14 74C 75C 7 1 2 8 76C (I) 77C 9 3 4 10 78C 79C 15 11 12 16 80C 81C The 12-point interpolation is not possible if either of the top 82C two rows is above the original field northern latitude. The 83C nearest neighbour is used if both rows are above, and a 4-pt 84C bilinear interpolation is used if the top row is above. 85C Similarily, if either of the bottom two rows is below the original 86C field southern latitude. 87C 88C 89C Externals 90C --------- 91C 92C INTLOG - Logs messages 93C JGETGG - Reads the definition of a gaussian grid 94C HGENGRD - Calculates original lat/long (before rotation) for 95C a rotated grid 96C HNEI12 - Finds neighbours for points for interpolation 97C HWTS12 - Calculates weightings for points for interpolation 98C CHKPREC - Check if precipitation threshold has been redefined 99C FORCED_NEAREST_NEIGHBOUR - check forced interpolation method 100C HIRLAM_USERSPACE_1_GET - HIRLAM userspace memory #1 101C HIRLAM_USERSPACE_2_GET - HIRLAM userspace memory #2 102C HIRLAM_USERSPACE_3_GET - HIRLAM userspace memory #3 103C 104C 105C Reference 106C --------- 107C 108C None. 109C 110C 111C Comments 112C -------- 113C 114C None. 115C 116C 117C Author 118C ------ 119C 120C J.D.Chambers ECMWF January 2001 121C 122C 123C Modifications 124C ------------- 125C 126C None. 127C 128C----< 129C -----------------------------------------------------------------| 130C* Section 0. Definition of variables. 131C -----------------------------------------------------------------| 132C 133 IMPLICIT NONE 134C 135#include "parim.h" 136#include "nifld.common" 137#include "nofld.common" 138#include "grfixed.h" 139C 140C Parameters 141C 142 INTEGER JNORTH, JSOUTH, JWEST, JEAST, JW_E, JN_S 143 INTEGER JP12PT, JP4PT 144 PARAMETER (JP12PT = 0) 145 PARAMETER (JP4PT = 1) 146 PARAMETER (JNORTH = 1) 147 PARAMETER (JWEST = 2) 148 PARAMETER (JSOUTH = 3) 149 PARAMETER (JEAST = 4) 150 PARAMETER (JW_E = 1) 151 PARAMETER (JN_S = 2) 152C 153C Function arguments 154C 155 LOGICAL L12PNT 156 INTEGER KOUNT, KGAUSS, KSIZE, NLON, NLAT 157 REAL AREA(4), POLE(2), GRID(2) 158 REAL OLDFLD(*), NEWFLD(KSIZE) 159 CHARACTER*1 HTYPE 160C 161C Local variables 162C 163 INTEGER NEXT, LOOP, IRET, NLEN, NBYTES, NUMBER 164 INTEGER NEAREST, NEND, COUNT 165 INTEGER NJPR, NJPB 166 INTEGER NGSPEC_TMP 167 INTEGER IOGAUSS,IOREPR 168 CHARACTER*1 HTYPE_TMP 169C 170 INTEGER KSCHEME(1),NEIGH(12,1), KLA(1) 171 REAL PWTS(12,1) 172 POINTER (IPKSCHE, KSCHEME) 173 POINTER (IPNEIGH, NEIGH) 174 POINTER (IPKLA, KLA) 175 POINTER (IPPWTS, PWTS) 176 LOGICAL LVEGGY 177C 178 REAL PDLO0(1),PDLO1(1),PDLO2(1),PDLO3(1),PDLAT(1) 179 POINTER (IPPDLO0, PDLO0) 180 POINTER (IPPDLO1, PDLO1) 181 POINTER (IPPDLO2, PDLO2) 182 POINTER (IPPDLO3, PDLO3) 183 POINTER (IPPDLAT, PDLAT) 184C 185 INTEGER IGG 186 INTEGER KPTS(1) 187 REAL GLATS(1) 188 INTEGER IOFFS(1) 189 POINTER (IPKPTS, KPTS) 190 POINTER (IPIOFFS, IOFFS) 191 POINTER (IPGLATS, GLATS) 192C 193 INTEGER ILL 194 REAL RLAT(1),RLON(1) 195 POINTER (IPRLAT, RLAT) 196 POINTER (IPRLON, RLON) 197C 198 REAL OLD(KOUNT) 199C 200 CHARACTER*7 FNAME 201 PARAMETER (FNAME = 'HIRLAM ') 202C 203 DATA NUMBER/-1/,IOGAUSS/-1/,IOREPR/-1/ 204C 205 SAVE IPKSCHE, IPNEIGH, IPKLA, IPPWTS 206 SAVE IPPDLO0, IPPDLO1, IPPDLO2, IPPDLO3, IPPDLAT 207 SAVE IPKPTS, IPIOFFS, IPGLATS 208 SAVE IPRLAT, IPRLON 209 SAVE NUMBER, IOGAUSS, IOREPR 210C 211C Externals 212C 213 LOGICAL FORCED_NEAREST_NEIGHBOUR 214 LOGICAL ISOCTAHEDRAL 215 INTEGER HNEI12 216 INTEGER HGENGRD 217#ifdef POINTER_64 218 INTEGER*8 219#else 220 INTEGER 221#endif 222 X HIRLAM_USERSPACE_1_GET, 223 X HIRLAM_USERSPACE_2_GET, 224 X HIRLAM_USERSPACE_3_GET 225C 226C Statement functions 227C 228 REAL A, B 229 LOGICAL NOTEQ 230C 231C Tests that A is not (almost) equal to B 232C 233 NOTEQ(A,B) = (ABS((A)-(B)).GT.(ABS(A)*1E-3)) 234C 235C -----------------------------------------------------------------| 236C Section 1. Initialise. 237C -----------------------------------------------------------------| 238C 239 100 CONTINUE 240C 241 HIRLAM = 0 242C 243 CALL JDEBUG() 244C 245 IF( L12PNT ) THEN 246 CALL INTLOG(JP_DEBUG,FNAME//': 12-pt interpolation',JPQUIET) 247 ELSE 248 CALL INTLOG(JP_DEBUG,FNAME//': 4-pt interpolation',JPQUIET) 249 ENDIF 250C 251 CALL CHKPREC() 252 IF( LPREC )THEN 253 CALL INTLOG(JP_DEBUG, 254 X FNAME//': precipitation threshold applied',JPQUIET) 255 ELSE 256 CALL INTLOG(JP_DEBUG, 257 X FNAME//': precipitation threshold not applied',JPQUIET) 258 ENDIF 259 260 IF( HTYPE.NE.'O' .AND. HTYPE.NE.'R' .AND. HTYPE.NE.'U' 261 X .AND. HTYPE.NE.'F' ) THEN 262 CALL INTLOG(JP_ERROR, 263 X FNAME//': Gaussian type unsupported ('//HTYPE//')',JPQUIET) 264 HIRLAM = 1 265 GOTO 900 266 ENDIF 267 268C Set interpolation method (use nearest neighbour if required) 269 LVEGGY = FORCED_NEAREST_NEIGHBOUR(LMETHOD,NITABLE,NIPARAM) 270 IF( LVEGGY ) CALL INTLOG(JP_DEBUG, 271 X FNAME//': nearest neighbour processing',JPQUIET) 272 273 274C 275C Dynamically allocate memory for gaussian grid information. 276C 277 IGG = KGAUSS*2 278 NBYTES = (IGG*JPRLEN) + (2*IGG+1)*JPBYTES 279C 280 IPKPTS = HIRLAM_USERSPACE_1_GET(NBYTES) 281#ifdef hpR64 282 IPKPTS = IPKPTS/(1024*1024*1024*4) 283#endif 284 IF( IPKPTS.EQ.0 ) THEN 285 CALL INTLOG(JP_ERROR, 286 X FNAME//': Memory allocate fail',JPQUIET) 287 HIRLAM = 1 288 GOTO 900 289 ENDIF 290C 291 IPGLATS = IPKPTS + (IGG*JPBYTES) 292 IPIOFFS = IPGLATS + (IGG*JPRLEN) 293C 294 IOGAUSS = -1 295 IOREPR = -1 296C 297C Build up offsets to start of each latitude in original field(s) 298C 299 IF( (KGAUSS.NE.IOGAUSS).OR.(IOREPR.NE.NIREPR) ) THEN 300 CALL INTLOG(JP_DEBUG,FNAME//': trying Gaussian '//HTYPE,KGAUSS) 301 CALL JGETGG(KGAUSS,HTYPE,GLATS,KPTS,IRET) 302 IF( IRET.NE.0 ) THEN 303 CALL INTLOG(JP_ERROR, 304 X FNAME//': JGETGG failed to get gaussian data',JPQUIET) 305 HIRLAM = 2 306 GOTO 900 307 ENDIF 308C 309 IOFFS(1) = 1 310 DO LOOP = 2, (KGAUSS*2+1) 311 IOFFS(LOOP) = IOFFS(LOOP-1) + KPTS(LOOP-1) 312 ENDDO 313C 314C 315C Check whether number of gaussian points agrees with grid specification. 316C If not and grid is N80/N160 non-octahedral, try the 12-pt specification. 317C 318 NUMBER = (IOFFS(KGAUSS*2+1) - 1) 319 IF( NUMBER.NE.KOUNT ) THEN 320 IF( (KGAUSS.EQ.80).OR.(KGAUSS.EQ.160) .AND. 321 X (NIREPR.EQ.JPQUASI) .AND. 322 X (.NOT. ISOCTAHEDRAL(NIGAUSS,MILLEN)) ) THEN 323 NGSPEC_TMP = NGSPEC 324 NGSPEC = 12 325 HTYPE_TMP = 'F' 326 IF( NIREPR.EQ.JPQUASI ) HTYPE_TMP = 'R' 327 CALL INTLOG(JP_DEBUG, 328 X FNAME//': trying 12-pt Gaussian '//HTYPE_TMP,KGAUSS) 329 CALL JGETGG(KGAUSS,HTYPE_TMP,GLATS,KPTS,IRET) 330 NGSPEC = NGSPEC_TMP 331C 332 IF( IRET.NE.0 ) THEN 333 CALL INTLOG(JP_ERROR, 334 X FNAME//': JGETGG failed to get gaussian data',JPQUIET) 335 HIRLAM = 2 336 GOTO 900 337 ENDIF 338C 339 IOFFS(1) = 1 340 DO LOOP = 2, (KGAUSS*2+1) 341 IOFFS(LOOP) = IOFFS(LOOP-1) + KPTS(LOOP-1) 342 ENDDO 343C 344 NUMBER = (IOFFS(KGAUSS*2+1) - 1) 345 IF( NUMBER.NE.KOUNT ) THEN 346 CALL INTLOG(JP_ERROR,FNAME//': Given no. of points =',KOUNT) 347 CALL INTLOG(JP_ERROR,FNAME//': Expected no. of pts =',NUMBER) 348 HIRLAM = 4 349 GOTO 900 350 ENDIF 351C 352 ELSE 353C 354 CALL INTLOG(JP_ERROR,FNAME//': Given no. of points =',KOUNT) 355 CALL INTLOG(JP_ERROR,FNAME//': Expected no. of pts =',NUMBER) 356 HIRLAM = 4 357 GOTO 900 358 ENDIF 359 ENDIF 360C 361 IOGAUSS = KGAUSS 362 IOREPR = NIREPR 363 ENDIF 364C 365C Preserve the input field 366C (in case OLDFLD and NEWFLD are the same arrays) 367C 368 NUMBER = (IOFFS(KGAUSS*2+1) - 1) 369 OLD(1:NUMBER) = OLDFLD(1:NUMBER) 370C 371C -----------------------------------------------------------------| 372C Section 2. Generate the lat/long points for the output grid 373C -----------------------------------------------------------------| 374C 375 200 CONTINUE 376C 377C (DJ Mar-05) 378 NLON = 1 + NINT((AREA(JEAST) - AREA(JWEST)) / GRID(JW_E)) 379 NLAT = 1 + NINT((AREA(JNORTH) - AREA(JSOUTH)) / GRID(JN_S)) 380C 381 NLEN = NLON * NLAT 382 383 NOWE = NLON 384 NONS = NLAT 385C 386C Check that given array is big enough for the new field. 387C 388 IF( NLEN.GT.KSIZE ) THEN 389 CALL INTLOG(JP_ERROR,FNAME//': Given array size = ',KSIZE) 390 CALL INTLOG(JP_ERROR,FNAME//': Required size = = ',NLEN) 391 HIRLAM = 5 392 GOTO 900 393 ENDIF 394C 395C Dynamically allocate memory for lat/long arrays. 396C 397 ILL = NLEN 398 NBYTES = 2*ILL*JPRLEN 399C 400 IPRLON = HIRLAM_USERSPACE_2_GET(NBYTES) 401#ifdef hpR64 402 IPRLON = IPRLON/(1024*1024*1024*4) 403#endif 404 IF( IPRLON.EQ.0 ) THEN 405 CALL INTLOG(JP_ERROR, 406 X FNAME//': Memory allocate fail',JPQUIET) 407 HIRLAM = 6 408 GOTO 900 409 ENDIF 410C 411 IPRLAT = IPRLON + (ILL*JPRLEN) 412C 413 IRET = HGENGRD(AREA,POLE,GRID,NLON,NLAT,RLAT,RLON) 414 IF( IRET.NE.0 ) THEN 415 CALL INTLOG(JP_ERROR, 416 X FNAME//': HGENGRD failed to get lat/lon grid data',JPQUIET) 417 HIRLAM = 7 418 GOTO 900 419 ENDIF 420C 421C -----------------------------------------------------------------| 422C Section 3. Find neighbours for each point for interpolation. 423C -----------------------------------------------------------------| 424C 425 300 CONTINUE 426C 427C Dynamically allocate memory for interpolation arrays. 428C 429 NJPR=17 430 NJPB=14 431 NBYTES = (NJPR*JPRLEN + NJPB*JPBYTES) * ILL 432C 433 IPPDLO0 = HIRLAM_USERSPACE_3_GET(NBYTES) 434#ifdef hpR64 435 IPPDLO0 = IPPDLO0/(1024*1024*1024*4) 436#endif 437 IF( IPPDLO0.EQ.0 ) THEN 438 CALL INTLOG(JP_ERROR, 439 X FNAME//': Memory allocate fail',JPQUIET) 440 HIRLAM = 8 441 GOTO 900 442 ENDIF 443C 444 IPPDLO1 = IPPDLO0 + (ILL*JPRLEN) 445 IPPDLO2 = IPPDLO1 + (ILL*JPRLEN) 446 IPPDLO3 = IPPDLO2 + (ILL*JPRLEN) 447 IPPDLAT = IPPDLO3 + (ILL*JPRLEN) 448 IPPWTS = IPPDLAT + (ILL*JPRLEN) 449 IPKSCHE = IPPWTS + (12*ILL*JPRLEN) 450 IPKLA = IPKSCHE + (ILL*JPBYTES) 451 IPNEIGH = IPKLA + (ILL*JPBYTES) 452C 453C Find neighbours. 454C 455 IRET = HNEI12(L12PNT,NLEN,RLAT,RLON,KGAUSS,KPTS,GLATS, 456 X KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH) 457 IF( IRET.NE.0 ) THEN 458 CALL INTLOG(JP_ERROR, 459 X FNAME//': HNEI12 failed to find neighbours',JPQUIET) 460 HIRLAM = 9 461 GOTO 900 462 ENDIF 463C 464C -----------------------------------------------------------------| 465C Section 4. Perform the 12-point horizontal interpolation. 466C -----------------------------------------------------------------| 467C 468 400 CONTINUE 469C 470C Setup the 12-point horizontal interpolation weights 471C 472 CALL HWTS12 473 X (NLEN,KSCHEME,KLA,PDLAT,GLATS,PDLO0,PDLO1,PDLO2,PDLO3,NEIGH, 474 X PWTS) 475C 476C Calculate the interpolated grid point values 477C 478 DO LOOP = 1, NLEN 479 IF( KSCHEME(LOOP).EQ.JP12PT ) THEN 480C 481C See if any of the neighbours are missing 482C 483 COUNT = 0 484 IF( NOTEQ(OLD(NEIGH( 1,LOOP)),RMISSGV) ) COUNT = COUNT + 1 485 IF( NOTEQ(OLD(NEIGH( 2,LOOP)),RMISSGV) ) COUNT = COUNT + 1 486 IF( NOTEQ(OLD(NEIGH( 3,LOOP)),RMISSGV) ) COUNT = COUNT + 1 487 IF( NOTEQ(OLD(NEIGH( 4,LOOP)),RMISSGV) ) COUNT = COUNT + 1 488 IF( NOTEQ(OLD(NEIGH( 5,LOOP)),RMISSGV) ) COUNT = COUNT + 1 489 IF( NOTEQ(OLD(NEIGH( 6,LOOP)),RMISSGV) ) COUNT = COUNT + 1 490 IF( NOTEQ(OLD(NEIGH( 7,LOOP)),RMISSGV) ) COUNT = COUNT + 1 491 IF( NOTEQ(OLD(NEIGH( 8,LOOP)),RMISSGV) ) COUNT = COUNT + 1 492 IF( NOTEQ(OLD(NEIGH( 9,LOOP)),RMISSGV) ) COUNT = COUNT + 1 493 IF( NOTEQ(OLD(NEIGH(10,LOOP)),RMISSGV) ) COUNT = COUNT + 1 494 IF( NOTEQ(OLD(NEIGH(11,LOOP)),RMISSGV) ) COUNT = COUNT + 1 495 IF( NOTEQ(OLD(NEIGH(12,LOOP)),RMISSGV) ) COUNT = COUNT + 1 496C 497C Interpolate using twelve neighbours if none are missing 498C 499 IF( LVEGGY) THEN 500 NEAREST = 1 501 IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2 502 IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3 503 IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4 504 IF( PWTS( 5,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 5 505 IF( PWTS( 6,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 6 506 IF( PWTS( 7,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 7 507 IF( PWTS( 8,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 8 508 IF( PWTS( 9,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 9 509 IF( PWTS(10,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =10 510 IF( PWTS(11,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =11 511 IF( PWTS(12,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =12 512 NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP)) 513 ELSE 514 IF( COUNT.EQ.12) THEN 515 NEWFLD(LOOP) = 516 X OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) + 517 X OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) + 518 X OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) + 519 X OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP) + 520 X OLD(NEIGH( 5,LOOP)) * PWTS( 5,LOOP) + 521 X OLD(NEIGH( 6,LOOP)) * PWTS( 6,LOOP) + 522 X OLD(NEIGH( 7,LOOP)) * PWTS( 7,LOOP) + 523 X OLD(NEIGH( 8,LOOP)) * PWTS( 8,LOOP) + 524 X OLD(NEIGH( 9,LOOP)) * PWTS( 9,LOOP) + 525 X OLD(NEIGH(10,LOOP)) * PWTS(10,LOOP) + 526 X OLD(NEIGH(11,LOOP)) * PWTS(11,LOOP) + 527 X OLD(NEIGH(12,LOOP)) * PWTS(12,LOOP) 528C 529C Set missing if all neighbours are missing 530C 531 ELSE IF( COUNT.EQ.0 ) THEN 532 NEWFLD(LOOP) = RMISSGV 533C 534C Otherwise, use the nearest neighbour 535C 536 ELSE 537 NEAREST = 1 538 IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2 539 IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3 540 IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4 541 IF( PWTS( 5,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 5 542 IF( PWTS( 6,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 6 543 IF( PWTS( 7,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 7 544 IF( PWTS( 8,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 8 545 IF( PWTS( 9,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 9 546 IF( PWTS(10,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =10 547 IF( PWTS(11,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =11 548 IF( PWTS(12,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =12 549 NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP)) 550 ENDIF 551 ENDIF 552C 553 ELSE IF( KSCHEME(LOOP).EQ.JP4PT ) THEN 554C 555C See if any of the neighbours are missing 556C 557 COUNT = 0 558 IF( NOTEQ(OLD(NEIGH( 1,LOOP)),RMISSGV) ) COUNT = COUNT + 1 559 IF( NOTEQ(OLD(NEIGH( 2,LOOP)),RMISSGV) ) COUNT = COUNT + 1 560 IF( NOTEQ(OLD(NEIGH( 3,LOOP)),RMISSGV) ) COUNT = COUNT + 1 561 IF( NOTEQ(OLD(NEIGH( 4,LOOP)),RMISSGV) ) COUNT = COUNT + 1 562C 563C Interpolate using four neighbours if none are missing 564C 565c IF( COUNT.EQ.4.AND.LVEGGY.EQV.(.FALSE.) ) THEN 566 IF( LVEGGY) THEN 567 568 NEAREST = 1 569 IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2 570 IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3 571 IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4 572 NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP)) 573 574 ELSE 575 IF( COUNT.EQ.4) THEN 576 NEWFLD(LOOP) = 577 X OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) + 578 X OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) + 579 X OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) + 580 X OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP) 581C 582C Set missing if all neighbours are missing 583C 584 ELSE IF( COUNT.EQ.0 ) THEN 585 NEWFLD(LOOP) = RMISSGV 586C 587C Otherwise, use the nearest neighbour 588C 589 ELSE 590 NEAREST = 1 591 IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2 592 IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3 593 IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4 594 NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP)) 595 ENDIF 596 ENDIF 597C 598 ELSE 599 DO NEXT = 1, 4 600 IF( NEIGH(NEXT,LOOP).NE.0 ) 601 X NEWFLD(LOOP) = OLD(NEIGH(NEXT,LOOP)) 602 ENDDO 603C 604 ENDIF 605C 606C Remove precipitation if less than a 'trace' or if nearest 607C neighbour is less than a trace 608C 609 IF( LPREC ) THEN 610 IF( NEWFLD(LOOP).LT.ZPRECIP ) THEN 611 NEWFLD(LOOP) = 0.0 612 ELSE 613 NEAREST = 1 614 NEND = 12 615 IF( KSCHEME(LOOP).NE.JP12PT ) NEND = 4 616 DO NEXT = 2, NEND 617 IF( PWTS( NEXT,LOOP).GT.PWTS( NEAREST,LOOP) ) 618 X NEAREST = NEXT 619 ENDDO 620 IF( OLD(NEIGH(NEAREST,LOOP)).LT.ZPRECIP ) NEWFLD(LOOP) = 0.0 621 ENDIF 622 ENDIF 623C 624 ENDDO 625C 626C -----------------------------------------------------------------| 627C Section 9. Return. 628C -----------------------------------------------------------------| 629C 630 900 CONTINUE 631 RETURN 632 END 633 634