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 WVQLIN2(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE, 13 X NEWWAVE,NORTH,WEST,KPARAM,PMISS) 14C 15C----> 16C*****WVQLIN2* 17C 18C PURPOSE 19C ------- 20C 21C To interpolate a wave field quasi-regular latitude 22C longitude grid to a regular latitude/longitude grid. 23C 24C 25C INTERFACE 26C --------- 27C 28C IRET = WVQLIN2(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE, 29C X NEWWAVE,NORTH,WEST,KPARAM,PMISS) 30C 31C Input arguments 32C --------------- 33C 34C KNUM - No. of meridians from North to South pole (input field) 35C NUMPTS - Array giving number of points along each latitude 36C (empty latitudes have entry 0) 37C KE_W - First dimension of new array 38C = Number of points E-W in new grid 39C KN_S - Second dimension of new array 40C = Number of points N-S in new grid 41C RESON - Output grid resolution (degrees) 42C OLDWAVE - Original wave field 43C NORTH - Input and output grid northernmost latitude (degree) 44C WEST - Input and output grid westernmost longitude (degree) 45C KPARAM - Field parameter code 46C PMISS - Missing value indicator 47C 48C Output arguments 49C ---------------- 50C 51C NEWWAVE - New wave field 52C 53C Function returns 0 if the interpolation was OK. 54C 55C 56C METHOD 57C ------ 58C 59C A bi-linear interpolation is only done if all four neighbouring 60C points are sea. 61C The nearest neighbouring grid point value is taken if one of the 62C neighbouring points is missing. 63C The interpolated point value is set to the 'missing' indicator 64C value if the point is outside the input grid. 65C 66C A field representing a 'direction' is resolved into components 67C each of which is interpolated separately. 68C 69C The output area is the same as the input area. 70C 71C 72C EXTERNALS 73C --------- 74C 75C WVQLID2 - Calculate the indices of neighbouring points 76C INTLOG - Log error message. 77C 78C 79C REFERENCE 80C --------- 81C 82C Based on: 83C SUBROUTINE INTERPOLATE 84C Peter Janssen ECMWF September 1995 85C and: 86C SUBROUTINE EXPOINT 87C Heinz Gunther ECMWF December 1989 88C 89C 90C Author. 91C ------- 92C 93C J.D.Chambers ECMWF September 1996 94C 95C J.D.Chambers ECMWF September 1998 96C Modified to handle subarea input/output. 97C The output area is the same as the input area. 98C 99C J.D.Chambers ECMWF November 2000 100C Modified to input grid resolution to be different from the 101C output grid resolution. 102C 103C----< 104C 105 IMPLICIT NONE 106C 107#include "parim.h" 108C 109C Parameters 110C 111 INTEGER JPROUTINE, JPLLMAX, JPNMOUT, 112 X JPNW, JPNE, JPSW, JPSE, JPN, JPS, 113 X JPDISNW, JPDISNE, JPDISSW, JPDISSE 114 PARAMETER( JPROUTINE = 19400 ) 115 PARAMETER( JPLLMAX = 1801 ) ! allow up to 0.1 degree resolution 116 PARAMETER( JPNMOUT = 1800 ) ! allow up to 0.1 degree resolution 117 PARAMETER( JPNW=1, JPNE=2, JPSW=3, JPSE=4, JPN=5, JPS=6, 118 X JPDISNW=7, JPDISNE=8, JPDISSW=9, JPDISSE=10 ) 119C 120C Arguments 121C 122 INTEGER KNUM, NUMPTS(*), KE_W, KN_S, KPARAM 123 REAL RESON, NORTH, WEST, PMISS 124 REAL RNS !FIXME: Unused! Difference in degrees in NS disrection 125 REAL OLDWAVE(*), NEWWAVE(*) 126C 127C Local variables 128C 129 INTEGER NLAT, K, I, NEWROW, NEWCOL, I1N, I1S, I_N, I_S, I2N, I2S 130 REAL*4 DELONGN, DELONGS, DELAT, DI1N, DI1S, DI2N, DI2S 131 REAL*4 DIST_N, DIST_S 132 REAL*4 DISNW, DISNE, DISSW, DISSE 133 REAL*4 U1, U2, ZXIN, ZXIS 134 REAL LON(JPNMOUT*2), LAT(JPNMOUT+1) 135 LOGICAL INGRID, LDIREC 136 INTEGER NW_PT, NE_PT, SW_PT, SE_PT 137 REAL*4 NWVALUE, NEVALUE, SWVALUE, SEVALUE 138 REAL*4 RAD 139 REAL*4 CNW_PT, CNE_PT, CSW_PT, CSE_PT 140 REAL*4 SNW_PT, SNE_PT, SSW_PT, SSE_PT 141 REAL*4 C1, C2, S1, S2, CC, SS 142C 143 DATA RAD/0.017453293/ 144C 145 REAL OLDLATS(JPLLMAX), OLDEAST 146 REAL OLDNS 147 INTEGER ISTART, IEND 148 INTEGER LOCATE 149 INTEGER NEXT, NUMNEW(JPLLMAX) 150 INTEGER IRET 151 INTEGER XKNUM,XKE_W,XKN_S 152 DATA XKNUM/-1/,XKE_W/-1/,XKN_S/-1/ 153 REAL XRESON,XNORTH,XWEST 154 DATA XRESON/-999.0/,XNORTH/-999.0/,XWEST/-999.0/ 155 INTEGER IFIRST, NEWIDX(4,3600*1801) 156 DATA IFIRST/0/ 157 REAL*4 DISTNEW(10,3600*1801) 158C 159 INTEGER IOS 160 INTEGER OLDKNUM 161 DATA OLDKNUM/0/ 162 INTEGER INDEX(JPLLMAX) 163 REAL DEPS 164 DATA DEPS/0.00005/ 165C Inline function - tests TRUE if A is a missing value 166C EPS Tolerance on missing data flag 167 REAL EPS 168 DATA EPS/0.1/ 169 LOGICAL LSW,LSE,LNE,LNW,LTEMP 170C 171 SAVE IFIRST,XKNUM,XKE_W,XKN_S,XRESON,XNORTH,XWEST,NEWIDX,DISTNEW 172C 173C Externals 174C 175 INTEGER WVQLID2 176 LOGICAL IS_WAVE_DIRECTION 177 EXTERNAL WVQLID2, IS_WAVE_DIRECTION 178C 179C -----------------------------------------------------------------| 180C* Section 1. Initalisation. 181C -----------------------------------------------------------------| 182C 183 100 CONTINUE 184C 185 WVQLIN2 = 0 186 CALL INTLOG(JP_DEBUG, 187 X 'WVQLIN2: Wave interpolation requested.',JPQUIET) 188C 189C Check reduced latitude/longitude grid specification 190C 191 IF( KNUM.GT.JPLLMAX ) THEN 192 CALL INTLOG(JP_ERROR, 193 X 'WVQLIN2: Number of latitudes in input lat/long grid = ',KNUM) 194 CALL INTLOG(JP_ERROR, 195 X 'WVQLIN2: And is greater than allowed maximum = ',JPLLMAX) 196 WVQLIN2 = JPROUTINE + 1 197 GOTO 900 198 ENDIF 199C 200C Ensure working array dimensions are adequate for required output. 201C 202 IF( (KE_W.GT.JPNMOUT*2).OR.(KN_S.GT.JPNMOUT+1) ) THEN 203 CALL INTLOG(JP_ERROR, 204 X 'WVQLIN2: Internal array dimensions are too small',JPQUIET) 205 CALL INTLOG(JP_ERROR, 206 X 'WVQLIN2: for given lat/long output field.',JPQUIET) 207 WVQLIN2 = JPROUTINE + 3 208 GOTO 900 209 ENDIF 210C 211C Set up index for latitude lines in the input reduced lat/lon array 212C 213 INDEX(1) = 0 214 DO K = 2, KNUM 215 INDEX(K) = INDEX(K-1) + NUMPTS(K-1) 216 ENDDO 217 218 DO K = 1, KN_S 219 LAT(K) = NORTH - FLOAT(K-1)*RESON 220 NUMNEW(K) = KE_W 221 ENDDO 222C 223C Calculate latitudes and longitudes of output grid points. 224C 225 DO I = 1,KE_W 226 LON(I) = FLOAT(I-1)*RESON 227 ENDDO 228C 229C Calculate latitudes for input (old) grid. 230C (Input and output area are same) 231C 232 OLDNS = RNS/REAL(KNUM-1) 233 DO K = 1, KNUM 234 OLDLATS(K) = NORTH - FLOAT(K-1)*OLDNS 235 ENDDO 236C 237C Find first and last latitudes in input grid which are not empty 238C 239 ISTART = 1 240 DO K = 1, KNUM/2 241 IF( NUMPTS(K).EQ.0 ) ISTART = ISTART + 1 242 ENDDO 243C 244 IEND = KNUM 245 DO K = KNUM, (KNUM/2) + 1, -1 246 IF( NUMPTS(K).EQ.0 ) IEND = IEND - 1 247 ENDDO 248C 249C For the regular output grid, calculate its latitudes and fill 250C in the number of points along each latitude. 251C 252C 253C Setup East depending on whether or not there is wrap-around 254C 255 IF( ABS((KE_W*RESON)-360.0).LT.DEPS ) THEN 256 OLDEAST = WEST + KE_W*RESON 257 ELSE 258 OLDEAST = WEST + (KE_W-1)*RESON 259 ENDIF 260C 261C Initialise all points with 'missing data' indicator 262C 263 NEXT = 0 264 DO NEWROW = 1, KN_S 265 DO NEWCOL = 1, KE_W 266 NEXT = NEXT + 1 267 NEWWAVE(NEXT) = PMISS 268 ENDDO 269 ENDDO 270C 271C Wave direction parameters need special handling 272 LDIREC = IS_WAVE_DIRECTION(KPARAM) 273C 274C -----------------------------------------------------------------| 275C* Section 2. Interpolation. 276C -----------------------------------------------------------------| 277C 278 200 CONTINUE 279C 280C Only calculate the indices on the first time through 281C 282 IF( (IFIRST.EQ.0 ) .OR. 283 X ( KNUM.NE.XKNUM ) .OR. 284 X ( KE_W.NE.XKE_W ) .OR. 285 X ( KN_S.NE.XKN_S ) .OR. 286 X ( RESON.NE.XRESON ).OR. 287 X ( NORTH.NE.XNORTH ).OR. 288 X ( WEST.NE.XWEST ) ) THEN 289C 290 IRET = WVQLID2(KNUM,NUMPTS,OLDLATS,WEST,OLDEAST, 291 X KN_S,NUMNEW,LAT,NEWIDX,DISTNEW) 292 IFIRST = 1 293 XKNUM = KNUM 294 XKE_W = KE_W 295 XKN_S = KN_S 296 XRESON = RESON 297 XNORTH = NORTH 298 XWEST = WEST 299 ENDIF 300C 301 DELAT = OLDNS 302 NEXT = 0 303C 304 DO 220 NEWROW = 1, KN_S 305C 306C Set up the distance between new row and the old rows to N and S. 307C 308 NLAT = NINT((NORTH - LAT(NEWROW))/DELAT) + 1 309 I_N = NLAT 310 I_S = MIN(I_N+1,KNUM*2) 311 DIST_N = ((NORTH - DELAT*(I_N-1)) - LAT(NEWROW)) / DELAT 312 DIST_S = 1.0 - DIST_N 313C 314C Check if the new interpolated row lies between 2 old rows which 315C have data points 316C 317 IF( (NUMPTS(I_N).GT.0).AND.(NUMPTS(I_S).GT.0) ) THEN 318C 319C Yes, use the calculated indices 320C 321 DO 210 NEWCOL = 1, KE_W 322C 323 NEXT = NEXT + 1 324C 325 NW_PT = NEWIDX(JPNW,NEXT) 326 NE_PT = NEWIDX(JPNE,NEXT) 327 SW_PT = NEWIDX(JPSW,NEXT) 328 SE_PT = NEWIDX(JPSE,NEXT) 329 NWVALUE = OLDWAVE(NW_PT) 330 NEVALUE = OLDWAVE(NE_PT) 331 SWVALUE = OLDWAVE(SW_PT) 332 SEVALUE = OLDWAVE(SE_PT) 333C 334C Test if any one of the four neighbouring points is missing. 335C 336 LTEMP = .FALSE. 337 LSE = .FALSE. 338 IF(ABS(SEVALUE-PMISS) .LT. EPS) LSE = .TRUE. 339 LSW = .FALSE. 340 IF(ABS(SWVALUE-PMISS) .LT. EPS) LSW = .TRUE. 341 LNE = .FALSE. 342 IF(ABS(NEVALUE-PMISS) .LT. EPS) LNE = .TRUE. 343 LNW = .FALSE. 344 IF(ABS(NWVALUE-PMISS) .LT. EPS) LNW = .TRUE. 345C 346 IF( (NW_PT.EQ.0 .OR. LNW) .OR. 347 X (NE_PT.EQ.0 .OR. LNE) .OR. 348 X (SW_PT.EQ.0 .OR. LSW) .OR. 349 X (SE_PT.EQ.0 .OR. LSE) ) THEN 350C If so, take nearest grid point value. 351C 352 DISNW = DISTNEW(JPDISNW,NEXT) 353 DISNE = DISTNEW(JPDISNE,NEXT) 354 DISSW = DISTNEW(JPDISSW,NEXT) 355 DISSE = DISTNEW(JPDISSE,NEXT) 356C 357 IF( (DISNW.LE.DISNE).AND. 358 X (DISNW.LE.DISSW).AND. 359 X (DISNW.LE.DISSE) ) THEN 360cs NEWWAVE(NEXT) = NWVALUE 361 LOCATE = NW_PT 362 LTEMP = LNW 363C 364 ELSE IF( (DISNE.LE.DISSW).AND. 365 X (DISNE.LE.DISSE) ) THEN 366cs NEWWAVE(NEXT) = NEVALUE 367 LOCATE = NE_PT 368 LTEMP = LNE 369C 370 ELSE IF( (DISSW.LE.DISSE) ) THEN 371cs NEWWAVE(NEXT) = SWVALUE 372 LOCATE = SW_PT 373 LTEMP = LSW 374C 375 ELSE 376cs NEWWAVE(NEXT) = SEVALUE 377 LOCATE = SE_PT 378 LTEMP = LSE 379C 380 ENDIF 381 382 IF( .NOT.LTEMP.AND.LOCATE.NE.0) THEN 383c IF( LOCATE.NE.0) THEN 384 NEWWAVE(NEXT) = OLDWAVE(LOCATE) 385 ENDIF 386 387C 388 ELSE 389C 390C Use bi-linear interpolation from four 391C neighbouring sea points. 392C 393 DI1N = DISTNEW(JPNW,NEXT) 394 DI2N = DISTNEW(JPNE,NEXT) 395 DI1S = DISTNEW(JPSW,NEXT) 396 DI2S = DISTNEW(JPSE,NEXT) 397 DIST_N = DISTNEW(JPN,NEXT) 398 DIST_S = DISTNEW(JPS,NEXT) 399C 400 NW_PT = OLDWAVE(NW_PT) 401 NE_PT = OLDWAVE(NE_PT) 402 SW_PT = OLDWAVE(SW_PT) 403 SE_PT = OLDWAVE(SE_PT) 404C 405 IF( .NOT. LDIREC ) THEN 406 U1 = NWVALUE*DI2N + NEVALUE*DI1N 407 U2 = SWVALUE*DI2S + SEVALUE*DI1S 408 NEWWAVE(NEXT) = U1*DIST_S + U2*DIST_N 409C 410 ELSE 411C 412C Fields representing a 'direction': resolve into 413C components and interpolate each separately. 414C 415 CNW_PT = COS(NWVALUE*RAD) 416 CNE_PT = COS(NEVALUE*RAD) 417 CSW_PT = COS(SWVALUE*RAD) 418 CSE_PT = COS(SEVALUE*RAD) 419 SNW_PT = SIN(NWVALUE*RAD) 420 SNE_PT = SIN(NEVALUE*RAD) 421 SSW_PT = SIN(SWVALUE*RAD) 422 SSE_PT = SIN(SEVALUE*RAD) 423 C1 = CNW_PT*DI2N + CNE_PT*DI1N 424 C2 = CSW_PT*DI2S + CSE_PT*DI1S 425 CC = C1*DIST_S + C2*DIST_N 426 S1 = SNW_PT*DI2N + SNE_PT*DI1N 427 S2 = SSW_PT*DI2S + SSE_PT*DI1S 428 SS = S1*DIST_S + S2*DIST_N 429 IF( SS.LT.0.0 ) THEN 430 NEWWAVE(NEXT) = ATAN2(SS,CC)/RAD + 360.0 431 ELSE 432 NEWWAVE(NEXT) = ATAN2(SS,CC)/RAD 433 ENDIF 434 ENDIF 435 ENDIF 436 210 CONTINUE 437 ELSE 438 NEXT = NEXT + KE_W 439 ENDIF 440C 441 220 CONTINUE 442C 443C -----------------------------------------------------------------| 444C* Section 9. Closedown. 445C -----------------------------------------------------------------| 446C 447 900 CONTINUE 448C 449 IF( WVQLIN2.EQ.0 ) THEN 450 CALL INTLOG(JP_DEBUG, 451 X 'WVQLIN2: Wave interpolation successful.',JPQUIET) 452 ELSE 453 CALL INTLOG(JP_DEBUG, 454 X 'WVQLIN2: Wave interpolation failed.',JPQUIET) 455 ENDIF 456c call abort() 457C 458 RETURN 459 END 460