1C 2C Copyright 1981-2016 ECMWF. 3C 4C This software is licensed under the terms of the Apache Licence 5C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 6C 7C In applying this licence, ECMWF does not waive the privileges and immunities 8C granted to it by virtue of its status as an intergovernmental organisation 9C nor does it submit to any jurisdiction. 10C 11 12 INTEGER FUNCTION WAVEXX2(NPARAM, NUMLATS, NPTS, 13 X NLATS, STEPNS, STEPWE, 14 X OLDWAVE, NEWWAVE, NORTH, WEST, PMISS) 15C 16C----> 17C*****WAVEXX2* 18C 19C PURPOSE 20C ------- 21C 22C Interpolates wave fields (except 2D spectra). 23C 24C 25C INTERFACE 26C --------- 27C 28C IRET = WAVEXX2(NPARAM, NUMLATS, NPTS, 29C X NLATS, STEPNS, STEPWE, 30C X OLDWAVE, NEWWAVE, NORTH, WEST, PMISS) 31C 32C Input arguments 33C --------------- 34C 35C NPARAM - Parameter ID 36C NUMLATS - Input lat number north-south 37C NPTS - Array giving number of points along each latitude 38C for input field (empty latitudes have entry 0) 39C NLATS - Number of points N-S in new grid 40C STEPNS - Output grid north-south resolution (degrees) 41C STEPWE - Output grid west-east resolution (degrees) 42C OLDWAVE - Original wave field 43C NORTH - Output grid northernmost latitude (degrees) 44C WEST - Output grid westernmost longitude (degrees) 45C PMISS - Missing data value 46C 47C Output arguments 48C ---------------- 49C 50C NEWWAVE - New wave field 51C 52C Function returns 0 if the interpolation was OK. 53C 54C 55C METHOD 56C ------ 57C 58C Builds the index of neighbouring points for the output grid. 59C Then works through the output grid points, checking for subarea 60C boundaries and looking up neighbouring point values and weights 61C (which may be missing data). 62C 63C 64C EXTERNALS 65C --------- 66C 67C WAVEIDX - Determines which nearest neighbour points to use for 68C interpolating to new output grid point 69C NUMPTWE - Calculates number of grid points between west/east 70C area boundaries 71C INTLOG - Log error message 72C JMALLOC - Dynamically allocate memory (deallocation -JFREE- not used) 73C 74C 75C REFERENCE 76C --------- 77C 78C None. 79C 80C 81C Author. 82C ------- 83C 84C S. Curic ECMWF Jun 2009 85C 86C 87C 88C----< 89C 90 IMPLICIT NONE 91C 92C Parameters 93C 94#include "parim.h" 95C 96 INTEGER JPROUTINE, JPMXLAT, 97 X JPNW, JPNE, JPSW, JPSE, JPN, JPS, 98 X JPDISNW, JPDISNE, JPDISSW, JPDISSE 99 PARAMETER( JPROUTINE = 40100 ) 100 PARAMETER( JPMXLAT = 1801 ) ! allow up to 0.1 degree resolution 101 PARAMETER( JPNW=1, JPNE=2, JPSW=3, JPSE=4, JPN=5, JPS=6, 102 X JPDISNW=7, JPDISNE=8, JPDISSW=9, JPDISSE=10 ) 103C 104#include "nifld.common" 105#include "nofld.common" 106#include "grspace.h" 107C 108C Arguments 109C 110 INTEGER NPARAM 111 INTEGER NUMLATS, NPTS(*), NLATS 112 REAL STEPNS, STEPWE, OLDWAVE(*), NEWWAVE(*), NORTH, WEST, PMISS 113C 114C Local variables 115C 116 INTEGER IEOFSET, INDEX, ISTART, IWEST, IWOFSET, KNEWNUM, 117 X KOLDNUM, LOOP, MISSLAT, NCOL, NEXT, 118 X NEXTWV, NROW, NUMNEW(JPMXLAT) 119 REAL AWEST, EAST, NEWLATS(JPMXLAT), OLDLATS(JPMXLAT), 120 X ONORTH, OSOUTH, PTLAT, PTLONG, RLATINC, ROWINC, 121 X SOUTH, OEAST, OWEST 122 INTEGER INE, INW, ISE, ISW 123 LOGICAL LDIREC 124 REAL C1, C2, CC, CNE_PT, CNW_PT, CSE_PT, CSW_PT, DI1N, DI1S, 125 X DI2N, DI2S, DISNE, DISNW, DISSE, DISSW, DK1, DK2, NE_PT, 126 X NW_PT, RAD, S1, S2, SE_PT, SNE_PT, SNW_PT, SS, SSE_PT, 127 X SSW_PT, SW_PT, U1, U2 128C 129 DATA RAD/0.01745329251/ 130C 131C Large, resolution-dependent arrays 132C 133 LOGICAL IS_ALLOCATED_LARGEARRAY 134 DATA IS_ALLOCATED_LARGEARRAY /.FALSE./ 135 SAVE IS_ALLOCATED_LARGEARRAY 136C 137 INTEGER NEWIDX(4,JPARRAYDIM_WAVE) 138#ifndef _CRAYFTN 139#ifdef POINTER_64 140 INTEGER*8 P_NEWIDX 141#endif 142#endif 143 POINTER(P_NEWIDX,NEWIDX) 144 SAVE P_NEWIDX 145C 146 REAL*4 DISTNEW(10,JPARRAYDIM_WAVE) 147#ifndef _CRAYFTN 148#ifdef POINTER_64 149 INTEGER*8 P_DISTNEW 150#endif 151#endif 152 POINTER(P_DISTNEW,DISTNEW) 153 SAVE P_DISTNEW 154C 155C Externals 156C 157 INTEGER WAVEIDX, NUMPTWE 158#ifdef POINTER_64 159 INTEGER*8 JMALLOC 160#else 161 INTEGER JMALLOC 162#endif 163 LOGICAL IS_WAVE_DIRECTION 164 EXTERNAL WAVEIDX, NUMPTWE, JMALLOC, IS_WAVE_DIRECTION 165C 166C -----------------------------------------------------------------| 167C* Section 1. Initalisation. 168C -----------------------------------------------------------------| 169C 170 100 CONTINUE 171C 172 WAVEXX2 = 0 173C 174C Initialise the bitmap value lookup function 175C 176 MISSLAT = 0 177 178 ONORTH = FLOAT(NIAREA(1))/PPMULT 179 OSOUTH = FLOAT(NIAREA(3))/PPMULT 180 RLATINC = FLOAT(NIGRID(2))/PPMULT 181C 182C Calculate number of latitudes if grid had been full from 183C North pole to South pole 184C 185 IF( NUMLATS.GT.JPMXLAT ) THEN 186 CALL INTLOG(JP_ERROR, 187 X 'WAVEXX2: Number of latitudes in input grid = ',NUMLATS) 188 CALL INTLOG(JP_ERROR, 189 X 'WAVEXX2: And is greater than allowed maximum = ',JPMXLAT) 190 WAVEXX2 = JPROUTINE + 1 191 GOTO 900 192 ENDIF 193C 194C 195C Fill an array with the number of points at each latitude for the 196C input field. 197C 198 IF(NIREPR.EQ.JPREDLL) THEN 199C 200C Input field is a reduced latitude/longitude grid 201C 202C .. but it may be 'pseudo-gaussian' in layout 203C (ie global, symmetric about the equator but no latitude 204C at the equator) 205C 206 IF( (ONORTH.NE.90.0).AND.(MOD(NUMLATS,2).EQ.0) ) THEN 207C 208 ONORTH = FLOAT(NIAREA(1)) 209 OSOUTH = FLOAT(NIAREA(3)) 210 RLATINC = FLOAT(NIGRID(2)) 211 212 ENDIF 213C 214 DO LOOP = 1, NUMLATS 215 OLDLATS(LOOP) = 90.0 - (LOOP-1)*RLATINC 216 ENDDO 217C 218 ELSE 219C 220C Input field is a regular latitude/longitude grid 221C 222 DO LOOP = 1, NUMLATS 223 OLDLATS(LOOP) = ONORTH - (LOOP-1)*RLATINC 224 ENDDO 225C 226 NPTS(1:NUMLATS) = NIWE 227C 228 ENDIF 229C 230C Allocate large, resolution-dependent arrays 231 IF (.NOT. IS_ALLOCATED_LARGEARRAY) THEN 232 P_NEWIDX = JMALLOC( 4*JPARRAYDIM_WAVE*JPBYTES) ! 4*len*bytes 233 P_DISTNEW = JMALLOC(10*JPARRAYDIM_WAVE*4) ! 10*len*bytes 234 IF( (P_NEWIDX.EQ.0) .OR. (P_DISTNEW.EQ.0) ) THEN 235 CALL INTLOG(JP_ERROR,'WAVEXX2: JMALLOC fail',JPQUIET) 236 WAVEXX2 = JPROUTINE + 1 237 GOTO 900 238 ENDIF 239 IS_ALLOCATED_LARGEARRAY = .TRUE. 240 END IF 241C 242C -----------------------------------------------------------------| 243C* Section 2. Setup number of points at each latitude for the 244C output latitude/longitude field. 245C -----------------------------------------------------------------| 246C 247 200 CONTINUE 248C 249 IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN) ) THEN 250C 251C Reduced (quasi-regular) gaussian output 252C 253 KNEWNUM = NOGAUSS*2 254 DO LOOP = 1, KNEWNUM 255 NUMNEW(LOOP) = NOLPTS(LOOP) 256 NEWLATS(LOOP) = ROGAUSS(LOOP) 257 ENDDO 258C 259 ELSE IF( NOREPR.EQ.JPREDLL ) THEN 260C 261C Reduced (quasi-regular) lat/long output 262C 263 KNEWNUM = NOREDLL 264 DO LOOP = 1, KNEWNUM 265 NUMNEW(LOOP) = NOLPTS(LOOP) 266 NEWLATS(LOOP) = ROREDLL(LOOP) 267 ENDDO 268C 269 ELSE 270C 271C Regular output 272C 273 MISSLAT = NINT((90.0 - NORTH)/STEPNS) 274 DO LOOP = 1, MISSLAT 275 NUMNEW(LOOP) = 0 276 ENDDO 277 DO LOOP = 1, NLATS 278 NUMNEW(LOOP+MISSLAT) = NINT(360.0/STEPWE) 279 ENDDO 280C 281 KNEWNUM = MISSLAT + NLATS 282 DO LOOP = 1, KNEWNUM 283 NEWLATS(LOOP) = 90.0 - (LOOP-1)*STEPNS 284 ENDDO 285C 286 ENDIF 287C 288C -----------------------------------------------------------------| 289C* Section 3. Get the input GRIB field characteristics. 290C -----------------------------------------------------------------| 291C 292 300 CONTINUE 293C 294C Calculate the indices of the input grid points to be used for 295C the output points 296C 297 OWEST = FLOAT(NIAREA(2))/PPMULT 298 OEAST = FLOAT(NIAREA(4))/PPMULT 299 300 WAVEXX2 = WAVEIDX(NUMLATS,NPTS,OLDLATS,OWEST,OEAST, 301 X KNEWNUM, NUMNEW, NEWLATS, NEWIDX, DISTNEW) 302 IF( WAVEXX2.NE.0 ) THEN 303 CALL INTLOG(JP_ERROR, 304 X 'WAVEXX2: Unable to calculate output grid indices',JPQUIET) 305 WAVEXX2 = JPROUTINE + 3 306 GOTO 900 307 ENDIF 308C 309C Wave direction parameters need special handling 310 LDIREC = IS_WAVE_DIRECTION(NPARAM) 311C 312C -----------------------------------------------------------------| 313C* Section 4. Work through the output subarea. 314C -----------------------------------------------------------------| 315C 316 400 CONTINUE 317C 318C Fill in the wave spectra values 319C 320 NEXT = 0 321 NEXTWV = 0 322C 323 SOUTH = NOAREA(3)/PPMULT 324 EAST = NOAREA(4)/PPMULT 325 ISTART = 0 326C 327C Work down through latitudes from north to south. 328C 329 DO NROW = 1, KNEWNUM 330C 331C If inside north-south (subarea) boundaries .. 332C 333 IF( (NOREPR.EQ.JPGAUSSIAN).OR.(NOREPR.EQ.JPQUASI) ) THEN 334 PTLAT = ROGAUSS(NROW) 335 ELSE 336 PTLAT = 90.0 - (NROW-1)*STEPNS 337 ENDIF 338C 339 IF( (PTLAT.LE.NORTH).AND.(ABS(PTLAT-SOUTH).GT.-0.0005) ) THEN 340C 341C Calculate number of points between west boundary of area and 342C Greenwich 343C 344 ROWINC = 360.0/NUMNEW(NROW) 345C 346 IWEST = INT(WEST/ROWINC) 347 IF (ABS(WEST - IWEST*ROWINC).GT.0.000001) THEN 348 IF (WEST.GT.0.0) IWEST = IWEST + 1 349 IF (WEST.LT.0.0) IWEST = IWEST - 1 350 ENDIF 351 AWEST = IWEST * ROWINC 352 IWOFSET = NUMPTWE(AWEST,0.0,ROWINC) 353 IEOFSET = NUMPTWE(AWEST,EAST,ROWINC) 354C 355C Work through subarea longitudes from west to east. 356C 357 DO NCOL = 1, NUMNEW(NROW) 358 PTLONG = AWEST + (NCOL-1)*ROWINC 359 NEXT = NUMPTWE(AWEST,PTLONG,ROWINC) 360 IF( (NEXT.LE.IEOFSET).AND.(NEXT.GE.0) ) THEN 361C 362C .. and inside west-east (subarea) boundaries 363C 364 NEXT = 1 + NEXT - IWOFSET 365 IF( NEXT.LE.0) NEXT = NEXT + NUMNEW(NROW) 366 NEXT = NEXT + ISTART 367 NEXTWV = NEXTWV + 1 368C 369 INW = NEWIDX(JPNW,NEXT) 370 INE = NEWIDX(JPNE,NEXT) 371 ISW = NEWIDX(JPSW,NEXT) 372 ISE = NEWIDX(JPSE,NEXT) 373C 374C Test if any of the four neighbouring points is missing. 375C 376 IF( (OLDWAVE(INW).EQ.PMISS) .OR. 377 X (OLDWAVE(ISW).EQ.PMISS) .OR. 378 X (OLDWAVE(INE).EQ.PMISS) .OR. 379 X (OLDWAVE(ISE).EQ.PMISS) ) THEN 380 ENDIF 381 382 IF( (INW.EQ.0) .OR. (OLDWAVE(INW).EQ.PMISS) .OR. 383 X (ISW.EQ.0) .OR. (OLDWAVE(ISW).EQ.PMISS) .OR. 384 X (INE.EQ.0) .OR. (OLDWAVE(INE).EQ.PMISS) .OR. 385 X (ISE.EQ.0) .OR. (OLDWAVE(ISE).EQ.PMISS) ) THEN 386 IF( (OLDWAVE(INW).EQ.PMISS) .AND. 387 X (OLDWAVE(ISW).EQ.PMISS) .AND. 388 X (OLDWAVE(INE).EQ.PMISS) .AND. 389 X (OLDWAVE(ISE).EQ.PMISS) ) THEN 390 ENDIF 391C 392C If so, take nearest grid point value. 393C 394 DISNW = DISTNEW(JPDISNW,NEXT) 395 DISNE = DISTNEW(JPDISNE,NEXT) 396 DISSW = DISTNEW(JPDISSW,NEXT) 397 DISSE = DISTNEW(JPDISSE,NEXT) 398C 399 IF( (DISNW.LE.DISNE).AND. 400 X (DISNW.LE.DISSW).AND. 401 X (DISNW.LE.DISSE)) THEN 402 INDEX = INW 403C 404 ELSE IF( (DISNE.LE.DISSW).AND. 405 X (DISNE.LE.DISSE) ) THEN 406 INDEX = INE 407C 408 ELSE IF( (DISSW.LE.DISSE) ) THEN 409 INDEX = ISW 410C 411 ELSE 412 INDEX = ISE 413 ENDIF 414C 415 IF(INDEX.EQ.0.OR.(OLDWAVE(INDEX).EQ.PMISS)) THEN 416C 417C Nearest point is missing 418C 419 NEWWAVE(NEXTWV) = PMISS 420C 421 ELSE 422 NEWWAVE(NEXTWV) = OLDWAVE(INDEX) 423 ENDIF 424C 425 ELSE 426C 427C Use bi-linear interpolation from four 428C neighbouring sea points (possible conv REAL*4 to REAL*8) 429C 430 DI1N = DISTNEW(JPNW,NEXT) 431 DI2N = DISTNEW(JPNE,NEXT) 432 DI1S = DISTNEW(JPSW,NEXT) 433 DI2S = DISTNEW(JPSE,NEXT) 434 DK1 = DISTNEW(JPN,NEXT) 435 DK2 = DISTNEW(JPS,NEXT) 436C 437 NW_PT = OLDWAVE(INW) 438 NE_PT = OLDWAVE(INE) 439 SW_PT = OLDWAVE(ISW) 440 SE_PT = OLDWAVE(ISE) 441 IF( .NOT. LDIREC ) THEN 442 U1 = NW_PT*DI2N + NE_PT*DI1N 443 U2 = SW_PT*DI2S + SE_PT*DI1S 444 NEWWAVE(NEXTWV) = U1*DK2 + U2*DK1 445 ELSE 446C 447C Fields representing a 'direction': resolve into 448C components and interpolate each separately. 449C 450#if (defined REAL_8) 451#define __ATAN2 DATAN2 452#define __COS DCOS 453#define __SIN DSIN 454#else 455#define __ATAN2 ATAN2 456#define __COS COS 457#define __SIN SIN 458#endif 459 CNW_PT = __COS(NW_PT*RAD) 460 CNE_PT = __COS(NE_PT*RAD) 461 CSW_PT = __COS(SW_PT*RAD) 462 CSE_PT = __COS(SE_PT*RAD) 463 SNW_PT = __SIN(NW_PT*RAD) 464 SNE_PT = __SIN(NE_PT*RAD) 465 SSW_PT = __SIN(SW_PT*RAD) 466 SSE_PT = __SIN(SE_PT*RAD) 467 C1 = CNW_PT*DI2N + CNE_PT*DI1N 468 C2 = CSW_PT*DI2S + CSE_PT*DI1S 469 CC = C1*DK2 + C2*DK1 470 S1 = SNW_PT*DI2N + SNE_PT*DI1N 471 S2 = SSW_PT*DI2S + SSE_PT*DI1S 472 SS = S1*DK2 + S2*DK1 473 IF( SS.LT.0.0 ) THEN 474 NEWWAVE(NEXTWV) = __ATAN2(SS,CC)/RAD + 360.0 475 ELSE 476 NEWWAVE(NEXTWV) = __ATAN2(SS,CC)/RAD 477 ENDIF 478#undef __SIN 479#undef __COS 480#undef __ATAN2 481 ENDIF 482 ENDIF 483 ENDIF 484 ENDDO 485C 486 ENDIF 487 ISTART = ISTART + NUMNEW(NROW) 488 ENDDO 489C 490C -----------------------------------------------------------------| 491C* Section 9. Closedown. 492C -----------------------------------------------------------------| 493C 494 900 CONTINUE 495 RETURN 496 END 497 498