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 PROGRAM TEST_INTERPOLATION 12 IMPLICIT NONE 13 14 15C Parameters 16 INTEGER JPGRIB ! GRIB size 17C PARAMETER (JPGRIB = 2000000) 18C PARAMETER (JPGRIB = 4284076) ! for R640/O640 output 19C PARAMETER (JPGRIB = 7000000) ! for INTUVP2 (requires 2 fields) 20 PARAMETER (JPGRIB = 33190420) ! for R1280/O1280 = 1/16 deg (INTUVP2 not tested) 21 INTEGER JPBYTES ! bytes/integer (8 on 64-bit architecture) 22#ifdef INTEGER_8 23 PARAMETER (JPBYTES = 8) 24#else 25 PARAMETER (JPBYTES = 4) 26#endif 27 LOGICAL IN, OUT 28 PARAMETER (IN=.TRUE., OUT=.FALSE.) 29 30 31C Locals 32 INTEGER R, I 33 INTEGER IO ! 0:?, 1:INTIN, 2: INTOUT 34 CHARACTER*256 ARG, KEY, VAL, F1, F2 35 CHARACTER*1 N 36 PARAMETER (N = CHAR(10)) 37 38 39C Internals & externals 40 LOGICAL SPLITARGS 41 INTEGER WRAP_INT, WRAP_INTF2, WRAP_INTUVP2, WRAP_INTVECT2 42 INTEGER IARGC 43 44 45C Initialize 46 R = 0 47 F1 = '' 48 F2 = '/dev/null' 49 50 51C Process command line, in order 52 DO I = 1, IARGC() 53 CALL GETARG(I,ARG) 54 R = -1 55 IO = 0 56 IF (SPLITARGS(ARG,IO,KEY,VAL)) THEN 57 58 IF (IO.EQ.0) THEN 59 IF (KEY.EQ.'input') F1 = VAL 60 IF (KEY.EQ.'output') F2 = VAL 61 IF (KEY.EQ.'input' .OR. KEY.EQ.'output') R=0 62 IF (KEY.EQ.'intf2') R=WRAP_INTF2 (F1,F2,JPGRIB,JPBYTES) 63 IF (KEY.EQ.'intuvp2') R=WRAP_INTUVP2 (F1,F2,JPGRIB,JPBYTES) 64 IF (KEY.EQ.'intvect2') R=WRAP_INTVECT2(F1,F2,JPGRIB,JPBYTES) 65 ELSEIF (IO.EQ.1) THEN 66 IF (KEY.EQ.'accuracy') R=WRAP_INT(IN, KEY,VAL,1,0) 67 IF (KEY.EQ.'area') R=WRAP_INT(IN, KEY,VAL,0,4) 68 IF (KEY.EQ.'date') R=WRAP_INT(IN, KEY,VAL,1,0) 69 IF (KEY.EQ.'form') R=WRAP_INT(IN, KEY,VAL,0,0) 70 IF (KEY.EQ.'grid') R=WRAP_INT(IN, KEY,VAL,0,2) 71 IF (KEY.EQ.'gridname') R=WRAP_INT(IN, KEY,VAL,0,0) 72 IF (KEY.EQ.'level') R=WRAP_INT(IN, KEY,VAL,1,0) 73 IF (KEY.EQ.'local') R=WRAP_INT(IN, KEY,VAL,1,0) 74 IF (KEY.EQ.'lsm_param') R=WRAP_INT(IN, KEY,VAL,0,0) 75 IF (KEY.EQ.'lsmset') R=WRAP_INT(IN, KEY,VAL,0,0) 76 IF (KEY.EQ.'matrix') R=WRAP_INT(IN, KEY,VAL,1,0) 77 IF (KEY.EQ.'npts') R=WRAP_INT(IN, KEY,VAL,2,0) 78 IF (KEY.EQ.'parameter') R=WRAP_INT(IN, KEY,VAL,1,0) 79 IF (KEY.EQ.'reduced') R=WRAP_INT(IN, KEY,VAL,1,0) 80 IF (KEY.EQ.'regular') R=WRAP_INT(IN, KEY,VAL,1,0) 81 IF (KEY.EQ.'scan') R=WRAP_INT(IN, KEY,VAL,1,0) 82 IF (KEY.EQ.'table') R=WRAP_INT(IN, KEY,VAL,1,0) 83 IF (KEY.EQ.'truncation') R=WRAP_INT(IN, KEY,VAL,1,0) 84 IF (KEY.EQ.'uselsm') R=WRAP_INT(IN, KEY,VAL,0,0) 85 IF (KEY.EQ.'useprecip') R=WRAP_INT(IN, KEY,VAL,0,0) 86 IF (KEY.EQ.'usewind') R=WRAP_INT(IN, KEY,VAL,0,0) 87 IF (KEY.EQ.'user_reduced_gaussian'.OR. 88 . KEY.EQ.'user_regular_gaussian') 89 . R=WRAP_INT(IN, KEY,VAL,1,0) 90 ELSEIF (IO.EQ.2) THEN 91 IF (KEY.EQ.'12-point') R=WRAP_INT(OUT,KEY,VAL,0,0) 92 IF (KEY.EQ.'accuracy') R=WRAP_INT(OUT,KEY,VAL,1,0) 93 IF (KEY.EQ.'area') R=WRAP_INT(OUT,KEY,VAL,0,4) 94 IF (KEY.EQ.'autoresol') R=WRAP_INT(OUT,KEY,VAL,1,0) 95 IF (KEY.EQ.'bitmap') R=WRAP_INT(OUT,KEY,VAL,0,0) 96 IF (KEY.EQ.'form') R=WRAP_INT(OUT,KEY,VAL,0,0) 97 IF (KEY.EQ.'frame') R=WRAP_INT(OUT,KEY,VAL,1,0) 98 IF (KEY.EQ.'gauss_from_file') R=WRAP_INT(OUT,KEY,VAL,0,0) 99 IF (KEY.EQ.'gaussian') R=WRAP_INT(OUT,KEY,VAL,1,0) 100 IF (KEY.EQ.'grid') R=WRAP_INT(OUT,KEY,VAL,0,2) 101 IF (KEY.EQ.'gridname') R=WRAP_INT(OUT,KEY,VAL,0,0) 102 IF (KEY.EQ.'interpolation') R=WRAP_INT(OUT,KEY,VAL,0,0) 103 IF (KEY.EQ.'level') R=WRAP_INT(OUT,KEY,VAL,1,0) 104 IF (KEY.EQ.'levtype') R=WRAP_INT(OUT,KEY,VAL,1,0) 105 IF (KEY.EQ.'packing') R=WRAP_INT(OUT,KEY,VAL,0,0) 106 IF (KEY.EQ.'parameter') R=WRAP_INT(OUT,KEY,VAL,1,0) 107 IF (KEY.EQ.'red_latlon') R=WRAP_INT(OUT,KEY,VAL,0,2) 108 IF (KEY.EQ.'reduced') R=WRAP_INT(OUT,KEY,VAL,1,0) 109 IF (KEY.EQ.'regular') R=WRAP_INT(OUT,KEY,VAL,1,0) 110 IF (KEY.EQ.'rotation') R=WRAP_INT(OUT,KEY,VAL,0,2) 111 IF (KEY.EQ.'scan') R=WRAP_INT(OUT,KEY,VAL,1,0) 112 IF (KEY.EQ.'specification') R=WRAP_INT(OUT,KEY,VAL,1,0) 113 IF (KEY.EQ.'style') R=WRAP_INT(OUT,KEY,VAL,0,0) 114 IF (KEY.EQ.'table') R=WRAP_INT(OUT,KEY,VAL,1,0) 115 IF (KEY.EQ.'truncation') R=WRAP_INT(OUT,KEY,VAL,1,0) 116 IF (KEY.EQ.'intermediate_gaussian'.OR. 117 . KEY.EQ.'user_reduced_gaussian'.OR. 118 . KEY.EQ.'user_reduced_latlon' .OR. 119 . KEY.EQ.'user_regular_gaussian') 120 . R=WRAP_INT(OUT,KEY,VAL,1,0) 121 ENDIF 122 123 ENDIF 124 IF (R.NE.0) THEN 125 PRINT *, 'Problem with option "'//ARG(1:INDEX(ARG,' ')-1)//'"' 126 PRINT *, ' key: "'//KEY(1:INDEX(KEY,' ')-1)//'"' 127 PRINT *, ' value: "'//VAL(1:INDEX(VAL,' ')-1)//'"' 128 CALL EXIT(1) 129 ENDIF 130 ENDDO 131 132 IF (F1(1:1).EQ.' ' .OR. F2(1:1).EQ.' ') THEN 133 666 PRINT *,N 134 .//'Usage: int --input=1.grib --output=2.grib [options] ...'//N 135 .//'where in/output files and interpolation fn. are mandatory:'//N 136 .//' --input=C input file'//N 137 .//' --output=C output file (default: /dev/null)'//N 138 .//' --intf2: general/MARS interpolation'//N 139 .//' (allows multiple GRIB input)'//N 140 .//' --intuvp2: u/v fields interpolation'//N 141 .//' (requires vo/d or u/v pair on GRIB input)'//N 142 .//' --intvect2: u/v fields interpolation'//N 143 .//' (requires u/v pair on GRIB input)'//N 144 .//N 145 .//'and any of the INTIN/OUT options, interpreted in order'//N 146 .//'(some combinations are inconsistent):'//N 147 .//' --INTIN:truncation=I: SH truncation (input)'//N 148 .//' --INTOUT:12-point=C: 12-point bilinear int. ([yes|no])'//N 149 .//' --INTOUT:accuracy=I: pck. accuracy ([16|24|...)'//N 150 .//' --INTOUT:area=R/R/R/R: geographical area (N/W/S/E)'//N 151 .//' --INTOUT:autoresol=I: SH auto truncation on/off ([1|0])'//N 152 .//' --INTOUT:grid=R/R: regular lat-lon, or'//N 153 .//' --INTOUT:gridname=[F|N|O]I: specific grid, with:'//N 154 .//' F: regular Gaussian'//N 155 .//' N: reduced Gaussian/"quasi-regular"'//N 156 .//' O: reduced Gaussian/"octahedral"'//N 157 .//' I: points between equator and poles (N)'//N 158 .//' --INTOUT:interpolation=C: interpol. method (ref. code)'//N 159 .//' --INTOUT:packing=C: pck. method ([second|simple|...])'//N 160 .//' --INTOUT:reduced=I: reduced Gaussian grid (alternate)'//N 161 .//' --INTOUT:regular=I: regular Gaussian grid (alternate)'//N 162 .//' --INTOUT:rotation=R/R: South pole new position (lat/lon)'//N 163 .//' --INTOUT:style=C: style ([dissemination|...])'//N 164 .//' --INTOUT:truncation=I: SH truncation (output)'//N 165 .//N 166 .//'with the following abbreviations:'//N 167 .//' I: integer number'//N 168 .//' R: real number'//N 169 .//' C: text string'//N 170 STOP 171 ENDIF 172 173 174 STOP 175 ENDPROGRAM 176 177 178C ------------------------------------------------------------------ 179 LOGICAL FUNCTION SPLITARGS(ARG,IN,KEY,VAL) 180 IMPLICIT NONE 181 CHARACTER*256 ARG, KEY, VAL 182 INTEGER IN 183 184 IN = 0 185 KEY = '' 186 VAL = '' 187 IF (ARG(1:2).EQ.'--') THEN 188 KEY = ARG(3:) 189 IF (INDEX(ARG,'=').GT.3) VAL = ARG(INDEX(ARG,'=')+1:) 190 IF (INDEX(ARG,'=').GT.3) KEY = ARG(3:INDEX(ARG,'=')-1) 191 ENDIF 192 IF (KEY(1:6).EQ.'INTIN:') THEN 193 IN = 1 194 KEY = KEY(7:) 195 ELSEIF (KEY(1:7).EQ.'INTOUT:') THEN 196 IN = 2 197 KEY = KEY(8:) 198 ENDIF 199 SPLITARGS = (INDEX(KEY,' ').GT.1) 200 ENDFUNCTION 201 202 203C ------------------------------------------------------------------ 204 INTEGER FUNCTION WRAP_INT(IN,K,C,NI,NR) 205 IMPLICIT NONE 206 LOGICAL IN 207 CHARACTER*256 K, C 208 INTEGER NI, NR 209 210 CHARACTER*256 M 211 INTEGER I(4) 212 REAL R(4) 213 214 INTEGER INTIN, INTOUT 215 EXTERNAL INTIN, INTOUT 216 217 I = 0 218 R = 0. 219 IF (IN) THEN 220 M = 'intin' 221 ELSE 222 M = 'intout' 223 ENDIF 224 M = M(1:INDEX(M,' ')-1)//'("'//K(1:INDEX(K,' ')-1)//'")' 225 IF (NI.LT.0 .OR. NI.GT.4) NI = 0 226 IF (NR.LT.0 .OR. NR.GT.4) NR = 0 227 IF (NI.GE.1) THEN 228 DO WHILE (INDEX(C,'/').GT.0) 229 C(INDEX(C,'/'):INDEX(C,'/')) = ' ' 230 ENDDO 231 READ (C,*) I(1:NI) 232 C = '' 233 PRINT *, M(1:INDEX(M,' ')-1), I(1:NI) 234 ELSEIF (NR.GE.1) THEN 235 DO WHILE (INDEX(C,'/').GT.0) 236 C(INDEX(C,'/'):INDEX(C,'/')) = ' ' 237 ENDDO 238 READ (C,*) R(1:NR) 239 C = '' 240 PRINT *, M(1:INDEX(M,' ')-1), R(1:NR) 241 ELSE 242 PRINT *, M(1:INDEX(M,' ')-1), ' "'//C(1:INDEX(C,' ')-1)//'"' 243 ENDIF 244 IF (IN) THEN 245 WRAP_INT = INTIN (K,I,R,C) 246 ELSE 247 WRAP_INT = INTOUT(K,I,R,C) 248 ENDIF 249 IF (WRAP_INT.NE.0) PRINT *, M(1:INDEX(M,' ')-1)//' failed.' 250 ENDFUNCTION 251 252 253C ------------------------------------------------------------------ 254 INTEGER FUNCTION WRAP_INTF2(INFILE,OUTFILE,LENGRIB,LENBYTES) 255 IMPLICIT NONE 256 CHARACTER*256 INFILE, OUTFILE 257 INTEGER LENGRIB, LENBYTES 258 259 INTEGER INGRIB(LENGRIB), NEWFLD(LENGRIB) 260 INTEGER U1, U2, N, LENINT, IREC, IRET 261 262 INTEGER INTF2 263 EXTERNAL INTF2 264 265C Initialize 266 INGRIB = 0 267 NEWFLD = 0 268 CALL EXCEPT() 269 270C Open input and output files 271 IRET = 0 272 CALL PBOPEN(U1, INFILE, 'r', IRET) 273 CALL CHECK(IRET.NE.0, 'PBOPEN failed (r).') 274 CALL PBOPEN(U2, OUTFILE, 'w', IRET) 275 CALL CHECK(IRET.NE.0, 'PBOPEN failed (w).') 276 277C Start of loop on input fields 278 N = 0 279 120 CONTINUE 280 N = N + 1 281 282C Read next field 283 CALL PBGRIB(U1, INGRIB, LENGRIB*LENBYTES, IREC, IRET) 284 IF (IRET.EQ.-1) THEN 285 IRET = 0 286 GOTO 190 287 ENDIF 288 CALL CHECK(IRET.NE.0, 'PBGRIB failed.') 289 290C Interpolate 291 PRINT *, 'Interpolate field #', N 292 LENINT = LENGRIB 293 IRET = INTF2(INGRIB,IREC,NEWFLD,LENINT) 294 CALL CHECK(IRET.NE.0, 'INTF2 failed') 295 296C Write the new field to file 297 IF (LENINT.GT.0) THEN 298 CALL PBWRITE(U2, NEWFLD, LENINT, IRET) 299 ELSE 300 PRINT *, 'Output same as input' 301 CALL PBWRITE(U2, INGRIB, IREC, IRET) 302 ENDIF 303 CALL CHECK(IRET.LT.0, 'PBWRITE failed') 304 IRET = 0 305 306C Loop back for next field 307 IF (IRET.EQ.0) GOTO 120 308 309C Closedown 310 190 CONTINUE 311 WRAP_INTF2 = IRET 312 CALL PBCLOSE(U1, IRET) 313 CALL PBCLOSE(U2, IRET) 314 PRINT *, 'All done.' 315 ENDFUNCTION 316 317 318C ------------------------------------------------------------------ 319 INTEGER FUNCTION WRAP_INTUVP2(INFILE,OUTFILE,LENGRIB,LENBYTES) 320 IMPLICIT NONE 321 CHARACTER*256 INFILE, OUTFILE 322 INTEGER LENGRIB, LENBYTES 323 324 INTEGER IVOGRIB(LENGRIB), IDVGRIB(LENGRIB) 325 INTEGER IUGRIB(LENGRIB), IVGRIB(LENGRIB) 326 INTEGER U1, U2, N, LENINT, IREC, IRET 327 328 INTEGER INTUVP2 329 EXTERNAL INTUVP2 330 331C Open input and output files 332 CALL PBOPEN(U1, INFILE, 'r', IRET) 333 CALL CHECK(IRET.NE.0, 'PBOPEN failed (r).') 334 CALL PBOPEN(U2, OUTFILE, 'w', IRET) 335 CALL CHECK(IRET.NE.0, 'PBOPEN failed (w).') 336 337C Start of loop on input vo/d fields 338 N = 0 339 220 CONTINUE 340 N = N + 1 341 342C Read next vo/d 343 CALL PBGRIB(U1, IVOGRIB, LENGRIB*LENBYTES, IREC, IRET) 344 IF (IRET.EQ.-1) THEN 345 IRET = 0 346 GOTO 290 347 ENDIF 348 CALL CHECK(IRET.NE.0, 'PBGRIB failed (vo).') 349 CALL PBGRIB(U1, IDVGRIB, LENGRIB*LENBYTES, IREC, IRET) 350 IF (IRET.EQ.-1) THEN 351 IRET = 0 352 GOTO 290 353 ENDIF 354 CALL CHECK(IRET.NE.0, 'PBGRIB failed (d).') 355 356C Interpolate 357 PRINT *, 'Interpolate vo/d > u/v #', N 358 LENINT = LENGRIB 359 IRET = INTUVP2(IVOGRIB, IDVGRIB, LENGRIB, IUGRIB, IVGRIB, LENINT) 360 CALL CHECK(IRET.NE.0, 'INTUVP2 failed.') 361 362C Write the new u/v fields to file 363 IF (LENINT.GT.0) THEN 364 CALL PBWRITE(U2, IUGRIB, LENINT, IRET) 365 CALL CHECK(IRET.LT.LENINT, 'PBWRITE failed.') 366 CALL PBWRITE(U2, IVGRIB, LENINT, IRET) 367 CALL CHECK(IRET.LT.LENINT, 'PBWRITE failed.') 368 IRET = 0 369 ENDIF 370 371C Loop back for next field 372 IF (IRET.EQ.0) GOTO 220 373 374C Closedown 375 290 CONTINUE 376 WRAP_INTUVP2 = IRET 377 CALL PBCLOSE(U1, IRET) 378 CALL PBCLOSE(U2, IRET) 379 PRINT *, 'All done.' 380 ENDFUNCTION 381 382 383C ------------------------------------------------------------------ 384 INTEGER FUNCTION WRAP_INTVECT2(INFILE,OUTFILE,LENGRIB,LENBYTES) 385 IMPLICIT NONE 386 CHARACTER*256 INFILE, OUTFILE 387 INTEGER LENGRIB, LENBYTES 388 389 INTEGER IIAGRIB(LENGRIB), IIBGRIB(LENGRIB) 390 INTEGER IOAGRIB(LENGRIB), IOBGRIB(LENGRIB) 391 INTEGER U1, U2, N, LENINT, IREC, IRET 392 393 INTEGER INTVECT2 394 EXTERNAL INTVECT2 395 396C Open input and output files 397 CALL PBOPEN(U1, INFILE, 'r', IRET) 398 CALL CHECK(IRET.NE.0, 'PBOPEN failed (r).') 399 CALL PBOPEN(U2, OUTFILE, 'w', IRET) 400 CALL CHECK(IRET.NE.0, 'PBOPEN failed (w).') 401 402C Start of loop on input u/v fields 403 N = 0 404 220 CONTINUE 405 N = N + 1 406 407C Read next u/v 408 CALL PBGRIB(U1, IIAGRIB, LENGRIB*LENBYTES, IREC, IRET) 409 IF (IRET.EQ.-1) THEN 410 IRET = 0 411 GOTO 290 412 ENDIF 413 CALL CHECK(IRET.NE.0, 'PBGRIB failed (u).') 414 CALL PBGRIB(U1, IIBGRIB, LENGRIB*LENBYTES, IREC, IRET) 415 IF (IRET.EQ.-1) THEN 416 IRET = 0 417 GOTO 290 418 ENDIF 419 CALL CHECK(IRET.NE.0, 'PBGRIB failed (v).') 420 421C Interpolate 422 PRINT *, 'Interpolate u/v #', N 423 LENINT = LENGRIB 424 IRET = INTVECT2(IIAGRIB,IIBGRIB,LENGRIB,IOAGRIB,IOBGRIB,LENINT) 425 CALL CHECK(IRET.NE.0, 'INTVECT2 failed.') 426 427C Write the new u/v fields to file 428 IF (LENINT.GT.0) THEN 429 CALL PBWRITE(U2, IOAGRIB, LENINT, IRET) 430 CALL CHECK(IRET.LT.LENINT, 'PBWRITE failed.') 431 CALL PBWRITE(U2, IOBGRIB, LENINT, IRET) 432 CALL CHECK(IRET.LT.LENINT, 'PBWRITE failed.') 433 IRET = 0 434 ENDIF 435 436C Loop back for next field 437 IF (IRET.EQ.0) GOTO 220 438 439C Closedown 440 290 CONTINUE 441 WRAP_INTVECT2 = IRET 442 CALL PBCLOSE(U1, IRET) 443 CALL PBCLOSE(U2, IRET) 444 PRINT *, 'All done.' 445 ENDFUNCTION 446 447 448C ------------------------------------------------------------------ 449 SUBROUTINE CHECK(OOPS,MSG) 450 IMPLICIT NONE 451 LOGICAL OOPS 452 CHARACTER MSG*(*) 453 IF (OOPS) THEN 454 PRINT *, MSG(1:INDEX(MSG,' ')-1) 455 CALL EXIT(3) 456 ENDIF 457 ENDSUBROUTINE 458 459