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 WV2DIDX(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST, 12 X KNEWNUM,NUMNEW,NEWLATS,NEWWAVE) 13C 14C----> 15C*****WV2DIDX* 16C 17C PURPOSE 18C ------- 19C 20C Determines which nearest-neighbour values of an input global wave 21C 2D-spectra grid field to use for an output global wave 2D-spectra 22C grid field. 23C 24C 25C INTERFACE 26C --------- 27C 28C IRET = WV2DIDX(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST, 29C X KNEWNUM,NUMNEW,NEWLATS,NEWWAVE) 30C 31C Input arguments 32C --------------- 33C 34C KOLDNUM - No. of meridians from North to South pole (input field) 35C NUMOLD - Array giving number of points along each latitude 36C (empty latitudes have entry 0) 37C OLDLATS - input field latitudes 38C OLDWEST - western longitude of the input field (degrees) 39C OLDEAST - eastern longitude of the input field (degrees) 40C 41C KNEWNUM - No. of meridians from North to South pole (output field) 42C NUMNEW - Array giving number of points along each latitude 43C (empty latitudes have entry 0) 44C NEWLATS - output field latitudes 45C 46C Output arguments 47C ---------------- 48C 49C NEWWAVE - Indices of points to use 50C 51C Function returns 0 if the interpolation was OK. 52C 53C 54C METHOD 55C ------ 56C 57C The index of the nearest neighbouring grid point value of 58C the input field is assigned to each point of the output 59C grid. 60C The index is zero if the output grid point is 'missing'. 61C 62C The input field can be regular or quasi-regular, and can be 63C global or a subarea. 64C The output field can be regular or quasi-regular. 65C 66C 67C EXTERNALS 68C --------- 69C 70C INTLOG - Log error message. 71C 72C 73C REFERENCE 74C --------- 75C 76C None 77C 78C 79C Author. 80C ------- 81C 82C J.D.Chambers ECMWF October 1997 83C 84C 85C----< 86C 87 IMPLICIT NONE 88C 89#include "parim.h" 90C 91C Parameters 92C 93 INTEGER JPROUTINE 94 PARAMETER ( JPROUTINE = 19410 ) 95 INTEGER JPLLMAX 96 PARAMETER ( JPLLMAX = 1801 ) 97C `--> allow upto 0.1 degree resolution 98 INTEGER JPNMOUT 99 PARAMETER ( JPNMOUT = 1800 ) 100C `--> allow upto 0.1 degree resolution 101C 102C Function arguments 103C 104 INTEGER KOLDNUM, NUMOLD, KNEWNUM, NUMNEW 105 DIMENSION NUMOLD(*), NUMNEW(*) 106 REAL OLDWEST, OLDEAST, OLDLATS, NEWLATS 107 DIMENSION OLDLATS(*), NEWLATS(*) 108 INTEGER NEWWAVE 109 DIMENSION NEWWAVE(*) 110C 111C Local arguments 112C 113 INTEGER NEWCOL, NEWROW, INCOL, LOOP, NEXT, NUMMAX 114 INTEGER I_NW, I_SW, I_N, I_S, I_NE, I_SE 115 REAL*8 DELONGN, DELONGS, DIFF, RESOL, DINC, DEPS 116 REAL*8 DIST_NW, DIST_SW, DIST_NE, DIST_SE, DIST_N, DIST_S 117 REAL*4 DISNW, DISNE, DISSW, DISSE 118 REAL*8 ZXIN, ZXIS, LON, OWEST, OEAST 119 DIMENSION LON(JPNMOUT*2) 120 LOGICAL LINGNS, LINGWE 121 INTEGER INDEXI 122 DIMENSION INDEXI(JPLLMAX) 123 LOGICAL LGLOBAL, LFULLG, LGWRAP 124 DATA DEPS/0.00005/ 125C 126C --------------------------------------------------------------------- 127C* Section 1. Initalisation. 128C --------------------------------------------------------------------- 129C 130 100 CONTINUE 131C 132 WV2DIDX = 0 133 CALL INTLOG(JP_DEBUG, 134 X 'WV2DIDX: Wave interpolation requested.',JPQUIET) 135C 136C Check latitude/longitude grid specification 137 IF( KOLDNUM.GT.JPLLMAX ) THEN 138 CALL INTLOG(JP_ERROR, 139 X 'WV2DIDX: Number of latitudes in input grid = ',KOLDNUM) 140 CALL INTLOG(JP_ERROR, 141 X 'WV2DIDX: And is greater than allowed maximum = ',JPLLMAX) 142 WV2DIDX = JPROUTINE + 1 143 GOTO 900 144 ENDIF 145C 146C Set up INDEXI for latitude lines in input lat/lon array 147C 148 149 INDEXI(1) = 0 150 DO LOOP = 2, KOLDNUM 151 INDEXI(LOOP) = INDEXI(LOOP-1) + NUMOLD(LOOP-1) 152 ENDDO 153C 154 OEAST = OLDEAST 155 OWEST = OLDWEST 156 157C Check whether the input field is global or a subarea 158C 159 NUMMAX = NUMOLD(1) 160 DO LOOP = 2, KOLDNUM 161 IF( NUMOLD(LOOP).GT.NUMMAX ) NUMMAX = NUMOLD(LOOP) 162 ENDDO 163 IF( NUMMAX.LE.0 ) THEN 164 CALL INTLOG(JP_ERROR, 165 X 'WV2DIDX: Input wave field has no points',JPQUIET) 166 WV2DIDX = JPROUTINE + 2 167 GOTO 900 168 ENDIF 169C 170 LGWRAP = ABS((OEAST-OWEST) - DBLE(360.0)).LT.DEPS 171 DINC = 360.0/DBLE(NUMMAX) 172 LFULLG = ABS((DMOD((OEAST-OWEST+DBLE(360.0)),DBLE(360.0)) 173 X +DINC)-DBLE(360.0)).LT.DEPS 174C 175 LGLOBAL = LFULLG.OR.LGWRAP 176C 177C --------------------------------------------------------------------- 178C* Section 2. Interpolation. 179C --------------------------------------------------------------------- 180C 181 200 CONTINUE 182C 183 NEXT = 1 184C 185 DO 220 NEWROW = 1, KNEWNUM 186C 187C Find old latitude rows to north and south of new latitude and 188C set up the distance between new row and the old rows to N and S. 189C 190 DO LOOP = 1, KOLDNUM -1 191 IF( (NEWLATS(NEWROW).LE.OLDLATS(LOOP)) .AND. 192 X (NEWLATS(NEWROW).GT.OLDLATS(LOOP+1)) ) THEN 193 I_N = LOOP 194 I_S = MIN(LOOP+1,KOLDNUM) 195 GOTO 205 196 ENDIF 197 ENDDO 198C 199 IF( NEWLATS(NEWROW).EQ.OLDLATS(KOLDNUM) ) THEN 200 I_N = KOLDNUM 201 I_S = KOLDNUM 202 ELSE 203 I_N = -1 204 I_S = -1 205 ENDIF 206C 207 205 CONTINUE 208 IF( I_N.NE.I_S) THEN 209 DIFF = OLDLATS(I_N) - OLDLATS(I_S) 210 DIST_N = (OLDLATS(I_N) - NEWLATS(NEWROW)) / DIFF 211 ELSE 212 DIST_N = 1.0 213 ENDIF 214 DIST_S = 1.0 - DIST_N 215C 216C Check if the new interpolated row lies between 2 old rows which 217C have data points 218C 219 LINGNS = .FALSE. 220 IF( I_N.GT.0 ) THEN 221 IF( (NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).GT.0) ) THEN 222C 223C Yes, so set up the grid increments. 224C 225C 226C Note that the grid increments are different depending on 227C whether or not the grid is global; global grids have an 228C extra interval at the end (eg between 358.5 and 0.0), while 229C subareas have outer limits defined by 'west' and 'east' and 230C a given number of points in between. 231C 232 IF( LGLOBAL ) THEN 233 DELONGN = 360.0/DBLE(NUMOLD(I_N)) 234 DELONGS = 360.0/DBLE(NUMOLD(I_S)) 235 ELSE 236 DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/ 237 X DBLE(NUMOLD(I_N)) 238 239 DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0)) 240 X +DINC) / DBLE(NUMOLD(I_N)) 241 DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/ 242 X DBLE(NUMOLD(I_S)) 243 DELONGS = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0)) 244 X +DINC) / DBLE(NUMOLD(I_S)) 245 ENDIF 246c 247 248 LINGNS = .TRUE. 249 ENDIF 250C 251C The equator is given special treatment so that the northern 252C and southern hemispheres for the parameter 250 can (later) 253C be stitched together. 254C 255C If the input field finishs at the equator and the output 256C field has a line at the equator, use the input equator for 257C interpolation. 258C 259 IF( NEWLATS(NEWROW).EQ.0.0 .AND. 260 x ((NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).EQ.0)) ) THEN 261C 262C Yes, so set up the grid increments. 263C 264 IF( LGLOBAL ) THEN 265 DELONGN = 360.0/DBLE(NUMOLD(I_N)) 266 ELSE 267 DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/ 268 X DBLE(NUMOLD(I_N)) 269 DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0)) 270 X +DINC) / DBLE(NUMOLD(I_N)) 271 ENDIF 272 DELONGS = DELONGN 273 DIST_N = 0.0 274 DIST_S = 1.0 275 I_S = I_N 276 LINGNS = .TRUE. 277 ENDIF 278 ENDIF 279C 280C Setup longitudes for the current row 281C 282 INCOL = NUMNEW(NEWROW) 283 IF( INCOL.NE.0 ) THEN 284 RESOL = 360.0/DBLE(INCOL) 285 DO LOOP = 1, INCOL 286 LON(LOOP) = DBLE(LOOP-1)*RESOL 287 ENDDO 288C 289 DO 210 NEWCOL = 1, INCOL 290C 291 LINGWE = (LON(NEWCOL).GT.OWEST).OR.(LON(NEWCOL).LT.OEAST) 292C 293 IF( LINGNS .AND. LINGWE ) THEN 294C 295C If point to be interpolated is inside the input field , 296C work out distances from new point to four neighbours. 297C 298 ZXIN = DMOD( (LON(NEWCOL) + (360.0-OWEST) + DBLE(720.0)), 299 X DBLE(360.0)) 300cs ZXIN = LON(NEWCOL) 301 ZXIS = ZXIN 302 ZXIN = ZXIN/DELONGN + 1.0 303 ZXIS = ZXIS/DELONGS + 1.0 304 I_NW = INT(ZXIN) 305 I_NE = I_NW + 1 306 I_SW = INT(ZXIS) 307 I_SE = I_SW + 1 308C 309C Allow wrap-around 310C 311C 312 I_NW = MOD(I_NW,NUMOLD(I_N)) 313 I_NE = MOD(I_NE,NUMOLD(I_N)) 314 I_SW = MOD(I_SW,NUMOLD(I_S)) 315 I_SE = MOD(I_SE,NUMOLD(I_S)) 316 IF( I_NW.EQ.0 ) I_NW = NUMOLD(I_N) 317 IF( I_SW.EQ.0 ) I_SW = NUMOLD(I_S) 318 IF( I_NE.EQ.0 ) I_NE = NUMOLD(I_N) 319 IF( I_SE.EQ.0 ) I_SE = NUMOLD(I_S) 320C 321C Calculate distance from interpolated point to its neighbours 322C 323 DIST_NW = ZXIN - REAL(I_NW) 324 DIST_NE = 1.0 - DIST_NW 325 DIST_SW = ZXIS - REAL(I_SW) 326 DIST_SE = 1.0 - DIST_SW 327C 328C Take nearest grid point value. 329C 330 DISNW = DIST_NW*DIST_NW + DIST_N*DIST_N 331 DISNE = DIST_NE*DIST_NE + DIST_N*DIST_N 332 DISSW = DIST_SW*DIST_SW + DIST_S*DIST_S 333 DISSE = DIST_SE*DIST_SE + DIST_S*DIST_S 334C 335 IF( (DISNW.LE.DISNE).AND.(DISNW.LE.DISSW).AND. 336 X (DISNW.LE.DISSE) ) THEN 337C 338C Use north-west neighbour 339C 340 NEWWAVE(NEXT) = INDEXI(I_N)+I_NW 341C 342 ELSE IF( (DISNE.LE.DISSW).AND.(DISNE.LE.DISSE) ) THEN 343C 344C Use north-east neighbour 345C 346 NEWWAVE(NEXT) = INDEXI(I_N)+I_NE 347C 348 ELSE IF( DISSW.LE.DISSE ) THEN 349C 350C Use south-west neighbour 351C 352 NEWWAVE(NEXT) = INDEXI(I_S)+I_SW 353 ELSE 354C 355C Use south-east neighbour 356C 357 NEWWAVE(NEXT) = INDEXI(I_S)+I_SE 358 ENDIF 359C 360C If the point to be interpolated is outside the input field, 361C set it to 'land'. 362C 363 ELSE 364 NEWWAVE(NEXT) = 0 365 ENDIF 366C 367 NEXT = NEXT + 1 368C 369 210 CONTINUE 370C 371 ENDIF 372C 373 220 CONTINUE 374C 375C --------------------------------------------------------------------- 376C* Section 9. Closedown. 377C --------------------------------------------------------------------- 378C 379 900 CONTINUE 380C 381 IF( WV2DIDX.EQ.0 ) THEN 382 CALL INTLOG(JP_DEBUG, 'WV2DIDX: successful.',JPQUIET) 383 ELSE 384 CALL INTLOG(JP_ERROR, 'WV2DIDX: failed.',JPQUIET) 385 ENDIF 386C 387 RETURN 388 END 389