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 HNEI12 12 X (L12PNT,KLEN,RLAT,RLON,KGAUSS,KPTS,PLATIN, 13 X KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH) 14C 15C----> 16C**** HNEI12 17C 18C Purpose 19C ------- 20C 21C This routine accepts a vector of points and finds the 12 22C neighbours for each point suitable for horizontal interpolation. 23C 24C 25C Interface 26C --------- 27C 28C IRET = HNEI12(L12PNT,KLEN,RLAT,RLON,KGAUSS,KPTS,PLATIN, 29C X KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH) 30C 31C 32C Input parameters 33C ---------------- 34C 35C L12PNT - Chooses between 12-point and 4-point interpolation 36C = .TRUE. for 12-point horizontal 37C = .FALSE. for 4-point 38C KLEN - Number of points along the vector 39C RLAT - List of latitudes for points. 40C RLON - List of longitudes for points. 41C KGAUSS - Original gaussian field number. 42C KPTS - List giving number of points in each latitude in original 43C gaussian field. 44C PLATIN - Original gaussian field latitudes (north and south). 45C 46C 47C Output parameters 48C ----------------- 49C 50C KSCHEME - Flag showing interpolation scheme to use for point 51C 0 = 12-point 52C 1 = 4-point bilinear 53C 2 = nearest neighbour 54C PDLAT - Meridian linear weight. 55C PDLO0 - Zonal linear weight for the latitude of point 5. 56C PDLO1 - Zonal linear weight for the latitude of point 1. 57C PDLO2 - Zonal linear weight for the latitude of point 3. 58C PDLO3 - Zonal linear weight for the latitude of point 11. 59C KLA - Latitude number in original field of latitude north of 60C each point in the vector. 61C NEIGH - List of indices in the original field of neighbouring 62C point values for each point in the vector. 63C 64C Returns 0 if function successful, non-zero otherwise. 65C 66C Common block usage 67C ------------------ 68C 69C None. 70C 71C 72C Method 73C ------ 74C 75C Numbering of the points (I is the interpolation point): 76C 77C 13 5 6 14 78C 79C 7 1 2 8 80C (I) 81C 9 3 4 10 82C 83C 15 11 12 16 84C 85C The 12-point interpolation is not possible if either of the top 86C two rows is above the original field northern latitude. The 87C nearest neighbour is used if both rows are above, and a 4-pt 88C bilinear interpolation is used if the top row is above. 89C Similarily, if either of the bottom two rows is below the original 90C field southern latitude. 91C 92C 93C Externals 94C --------- 95C 96C INTLOG - Log error message. 97C JNORSGG - Finds gaussian latitude north or south of given latitude 98C JDEBUG - Tests if debug output required. 99C JMALLOC - Dynamically allocate memory 100C JFREE - Free dynamically allocated memory 101C 102C 103C Reference 104C --------- 105C 106C ECMWF Meteorological Bulletin M1.6/7 107C IFS Documentation 108C Part VI: Technical and Computational Procedures (CY21R4) 109C March 2000 110C Section 2.3 111C 112C 113C Comments 114C -------- 115C 116C None. 117C 118C 119C Author 120C ------ 121C 122C J.D.Chambers ECMWF January 2001 123C 124C 125C Modifications 126C ------------- 127C 128C None. 129C 130C----< 131C -----------------------------------------------------------------| 132C* Section 0. Definition of variables. 133C -----------------------------------------------------------------| 134C 135 IMPLICIT NONE 136#include "parim.h" 137C 138C Parameters 139C 140 INTEGER JP12PT, JP4PT, JPNEARN 141 PARAMETER (JP12PT = 0) 142 PARAMETER (JP4PT = 1) 143 PARAMETER (JPNEARN = 2) 144C 145C Function arguments 146C 147 LOGICAL L12PNT 148 INTEGER KLEN, KGAUSS, KPTS(KGAUSS*2) 149 REAL RLAT(KLEN), RLON(KLEN), PLATIN(KGAUSS*2) 150C 151 INTEGER KLA(KLEN), KSCHEME(KLEN) 152 REAL PDLAT(KLEN),PDLO0(KLEN),PDLO1(KLEN),PDLO2(KLEN),PDLO3(KLEN) 153 INTEGER NEIGH(12,*) 154C 155C Local variables 156C 157 INTEGER NEXT, LOOP, NBYTES 158 INTEGER NLAT0, NLAT1, NLAT2, NLAT3, NLAT(12), NLONG(12) 159 REAL STEP0, STEP1, STEP2, STEP3 160 INTEGER NPREV, ISIZE, ISIZOLD, IOFFSET(1) 161 POINTER( IPOFF, IOFFSET ) 162C 163 DATA ISIZOLD/-1/, NPREV/-1/ 164 SAVE ISIZOLD,IPOFF,NPREV 165C 166C Externals 167C 168 INTEGER JNORSGG 169#ifdef POINTER_64 170 INTEGER*8 JMALLOC 171#else 172 INTEGER JMALLOC 173#endif 174C 175C -----------------------------------------------------------------| 176C Section 1. Initialise. 177C -----------------------------------------------------------------| 178C 179 100 CONTINUE 180C 181 HNEI12 = 0 182C 183 CALL JDEBUG() 184C 185C Get memory for offset array. 186C 187 ISIZE = KGAUSS*2 + 1 188 IF( ISIZE.GT.ISIZOLD ) THEN 189C 190 IF( ISIZOLD.GT.0 ) CALL JFREE(IPOFF) 191C 192 NBYTES = ISIZE * JPBYTES 193 IPOFF = JMALLOC(NBYTES) 194#ifdef hpR64 195 IPOFF = IPOFF/(1024*1024*1024*4) 196#endif 197 IF( IPOFF.EQ.0 ) THEN 198 CALL INTLOG(JP_WARN,'HNEI12: Memory allocate fail ',JPQUIET) 199 HNEI12 = 1 200 GOTO 900 201 ENDIF 202C 203 ISIZOLD = ISIZE 204 NPREV = -1 205C 206 ENDIF 207C 208C Build up offsets to start of each latitude in the original field. 209C 210 IF( KGAUSS.NE.NPREV ) THEN 211 IOFFSET(1) = 1 212 DO LOOP = 2, (KGAUSS*2+1) 213 IOFFSET(LOOP) = IOFFSET(LOOP-1) + KPTS(LOOP-1) 214 ENDDO 215 NPREV = KGAUSS 216 ENDIF 217C 218C -----------------------------------------------------------------| 219C Section 2. Loop through points along the vector. 220C -----------------------------------------------------------------| 221C 222 200 CONTINUE 223C 224 DO NEXT = 1,KLEN 225C 226C Clear current neighbour indices 227C 228 DO LOOP = 1, 12 229 NEIGH(LOOP,NEXT) = 0 230 ENDDO 231C 232C Find latitude numbers for neighbours. 233C 234C First find latitude numbers to north and south of the point 235C 236C -----------------------------------------------------------------| 237C Section 3. Nearest neighbour to be used. 238C -----------------------------------------------------------------| 239C 240 300 CONTINUE 241C 242 IF( RLAT(NEXT).GT.PLATIN(1) ) THEN 243 CALL INTLOG(JP_DEBUG, 244 X 'HNEI12: Nearest neighbour to south used for point ',NEXT) 245C 246C Above northern latitude of original field, so use nearest 247C neighbour (point 3 or 4) 248C 249 STEP2 = 360.0 / REAL(KPTS(1)) 250 NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2) 251 IF( NLONG( 3).GT.KPTS(1) ) NLONG( 3) = 1 252 NLONG( 4) = NLONG(3) + 1 253 IF( NLONG( 4).GT.KPTS(1) ) NLONG( 4) = 1 254C 255 IF( ABS((RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 )) .LE. 256 X ABS((RLON(NEXT) - REAL(NLONG(4)-1)*STEP2 )) ) THEN 257 NEIGH(3,NEXT) = NLONG(3) 258 ELSE 259 NEIGH(4,NEXT) = NLONG(4) 260 ENDIF 261 KSCHEME(NEXT) = JPNEARN 262 GOTO 600 263 ENDIF 264C 265 IF( RLAT(NEXT).LT.PLATIN(KGAUSS*2) ) THEN 266 CALL INTLOG(JP_DEBUG, 267 X 'HNEI12: Nearest neighbour to north used for point ',NEXT) 268C 269C Below southern latitude of original field, so use nearest 270C neighbour (point 1 or 2) 271C 272 STEP1 = 360.0 / REAL(KPTS(KGAUSS*2)) 273 NLONG( 1) = 1 + INT(RLON(NEXT)/STEP1) 274 IF( NLONG( 1).GT.KPTS(KGAUSS*2) ) NLONG( 1) = 1 275 NLONG( 2) = NLONG(1) + 1 276 IF( NLONG( 2).GT.KPTS(KGAUSS*2) ) NLONG( 2) = 1 277C 278 IF( ABS((RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 )) .LE. 279 X ABS((RLON(NEXT) - REAL(NLONG(2)-1)*STEP1 )) ) THEN 280 NEIGH(1,NEXT) = IOFFSET(KGAUSS*2) + NLONG(1) - 1 281 ELSE 282 NEIGH(2,NEXT) = IOFFSET(KGAUSS*2) + NLONG(2) - 1 283 ENDIF 284 KSCHEME(NEXT) = JPNEARN 285 GOTO 600 286 ENDIF 287C 288 NLAT1 = JNORSGG(RLAT(NEXT),PLATIN,KGAUSS,1) 289 NLAT2 = JNORSGG(RLAT(NEXT),PLATIN,KGAUSS,0) 290C 291 IF( NLAT2.EQ.NLAT1) NLAT2 = NLAT1 + 1 292C 293C Find outer latitude numbers to north and south of the point 294C 295 NLAT0 = NLAT1 - 1 296 NLAT3 = NLAT2 + 1 297C 298C -----------------------------------------------------------------| 299C Section 4. 4-point bilinear to be used. 300C -----------------------------------------------------------------| 301C 302 400 CONTINUE 303C 304 IF( (NLAT0.EQ.0).OR.(NLAT3.GT.KGAUSS*2).OR.(.NOT.L12PNT) ) THEN 305Cjdc CALL INTLOG(JP_DEBUG, 306Cjdc X 'HNEI12: 4-pt bilinear interpolation used for point ',NEXT) 307C 308C Between two northern latitudes or two southern latitudes of 309C original field, so use 4-point bilinear interpolation 310C (points 1, 2, 3 and 4) 311C 312C Points 1 and 2 313C 314 STEP1 = 360.0 / REAL(KPTS(NLAT1)) 315 NLONG(1) = 1 + INT(RLON(NEXT)/STEP1) 316 IF( NLONG(1).GT.KPTS(NLAT1) ) NLONG(1) = 1 317 NLONG(2) = NLONG(1) + 1 318 IF( NLONG(2).GT.KPTS(NLAT1) ) NLONG(2) = 1 319C 320 NEIGH(1,NEXT) = IOFFSET(NLAT1) + NLONG(1) - 1 321 NEIGH(2,NEXT) = IOFFSET(NLAT1) + NLONG(2) - 1 322C 323C Points 3 and 4 324C 325 STEP2 = 360.0 / REAL(KPTS(NLAT2)) 326 NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2) 327 IF( NLONG( 3).GT.KPTS(NLAT2) ) NLONG( 3) = 1 328 NLONG( 4) = NLONG(3) + 1 329 IF( NLONG( 4).GT.KPTS(NLAT2) ) NLONG( 4) = 1 330C 331 NEIGH(3,NEXT) = IOFFSET(NLAT2) + NLONG(3) - 1 332 NEIGH(4,NEXT) = IOFFSET(NLAT2) + NLONG(4) - 1 333C 334C Calculate zonal linear weights for the latitudes. 335C 336 PDLO1(NEXT) = ( RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 ) / STEP1 337 PDLO2(NEXT) = ( RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 ) / STEP2 338C 339 KLA(NEXT) = NLAT1 340 PDLAT(NEXT) = ( RLAT(NEXT) - PLATIN(NLAT1) ) / 341 X ( PLATIN(NLAT2) - PLATIN(NLAT1) ) 342C 343 KSCHEME(NEXT) = JP4PT 344 GOTO 600 345 ENDIF 346C 347C -----------------------------------------------------------------| 348C Section 5. 12-point interpolation to be used. 349C -----------------------------------------------------------------| 350C 351 500 CONTINUE 352C 353 KSCHEME(NEXT) = JP12PT 354C 355 KLA(NEXT) = NLAT1 356 PDLAT(NEXT) = ( RLAT(NEXT) - PLATIN(NLAT1) ) / 357 X ( PLATIN(NLAT2) - PLATIN(NLAT1) ) 358C 359C Northernmost 360C 361 NLAT(5) = NLAT0 362 NLAT(6) = NLAT0 363C 364C Northern 365C 366 NLAT(7) = NLAT1 367 NLAT(1) = NLAT1 368 NLAT(2) = NLAT1 369 NLAT(8) = NLAT1 370C 371C Southern 372C 373 NLAT( 9) = NLAT2 374 NLAT( 3) = NLAT2 375 NLAT( 4) = NLAT2 376 NLAT(10) = NLAT2 377C 378C Southernmost 379C 380 NLAT(11) = NLAT3 381 NLAT(12) = NLAT3 382C 383C -----------------------------------------------------------------| 384C Find longitude numbers for neighbours. 385C -----------------------------------------------------------------| 386C 387C Northernmost 388C 389 STEP0 = 360.0 / REAL(KPTS(NLAT0)) 390 NLONG(5) = 1 + INT(RLON(NEXT)/STEP0) 391 IF( NLONG(5).GT.KPTS(NLAT0) ) NLONG(5) = 1 392 NLONG(6) = NLONG(5) + 1 393 IF( NLONG(6).GT.KPTS(NLAT0) ) NLONG(6) = 1 394C 395C Northern 396C 397 STEP1 = 360.0 / REAL(KPTS(NLAT1)) 398 NLONG(1) = 1 + INT(RLON(NEXT)/STEP1) 399 IF( NLONG(1).GT.KPTS(NLAT1) ) NLONG(1) = 1 400 NLONG(7) = NLONG(1) - 1 401 IF( NLONG(7).LT.1 ) NLONG(7) = KPTS(NLAT1) 402 NLONG(2) = NLONG(1) + 1 403 IF( NLONG(2).GT.KPTS(NLAT1) ) NLONG(2) = 1 404 NLONG(8) = NLONG(2) + 1 405 IF( NLONG(8).GT.KPTS(NLAT1) ) NLONG(8) = 1 406C 407C Southern 408C 409 STEP2 = 360.0 / REAL(KPTS(NLAT2)) 410 NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2) 411 IF( NLONG( 3).GT.KPTS(NLAT2) ) NLONG( 3) = 1 412 NLONG( 9) = NLONG(3) - 1 413 IF( NLONG( 9).LT.1 ) NLONG( 9) = KPTS(NLAT2) 414 NLONG( 4) = NLONG(3) + 1 415 IF( NLONG( 4).GT.KPTS(NLAT2) ) NLONG( 4) = 1 416 NLONG(10) = NLONG(4) + 1 417 IF( NLONG(10).GT.KPTS(NLAT2) ) NLONG(10) = 1 418C 419C Southernmost 420C 421 STEP3 = 360.0 / REAL(KPTS(NLAT3)) 422 NLONG(11) = 1 + INT(RLON(NEXT)/STEP3) 423 IF( NLONG(11).GT.KPTS(NLAT3) ) NLONG(11) = 1 424 NLONG(12) = NLONG(11) + 1 425 IF( NLONG(12).GT.KPTS(NLAT3) ) NLONG(12) = 1 426C 427C -----------------------------------------------------------------| 428C Calculate zonal linear weights for the latitudes. 429C -----------------------------------------------------------------| 430C 431 PDLO0(NEXT) = ( RLON(NEXT) - REAL(NLONG(5)-1)*STEP0 ) / STEP0 432C 433 PDLO1(NEXT) = ( RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 ) / STEP1 434C 435 PDLO2(NEXT) = ( RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 ) / STEP2 436C 437 PDLO3(NEXT) = ( RLON(NEXT) - REAL(NLONG(11)-1)*STEP3 ) / STEP3 438C 439C -----------------------------------------------------------------| 440C Store indices of the neighbours. 441C -----------------------------------------------------------------| 442C 443 DO LOOP = 1, 12 444 NEIGH(LOOP,NEXT) = IOFFSET(NLAT(LOOP)) + NLONG(LOOP) - 1 445 ENDDO 446C 447C -----------------------------------------------------------------| 448C Section 6. End of loop along vector of points. 449C -----------------------------------------------------------------| 450C 451 600 CONTINUE 452C 453 ENDDO 454C 455C -----------------------------------------------------------------| 456C Section 9. Return 457C -----------------------------------------------------------------| 458C 459 900 CONTINUE 460C 461 RETURN 462 END 463