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 INTUVXH(PUVSH,KNVALS,ZNFLDO,KUGRIB,KVGRIB, 12 X OUTLENU,OUTLENV) 13C 14C----> 15C**** INTUVXH 16C 17C Purpose 18C ------- 19C 20C Interpolate U or V component spectral fields to grid point. 21C 22C 23C Interface 24C --------- 25C 26C IRET = INTUVXH(PUVSH,KNVALS,ZNFLDO,KUGRIB,KVGRIB,OUTLENU,OUTLENV) 27C 28C Input 29C ----- 30C 31C PUVSH - Spectral U/V values. (U first) 32C KNVALS - Number of values in each wind component field. 33C ZNFLDO - Work array. 34C 35C 36C Output 37C ------ 38C 39C KUGRIB - Output wind U component field (GRIB format). 40C KVGRIB - Output wind V component field (GRIB format). 41C OUTLENU - Output U field length (words). 42C OUTLENV - Output V field length (words). 43C 44C 45C Method 46C ------ 47C 48C None. 49C 50C 51C Externals 52C --------- 53C 54C INTUVDH - Encode/decode data into/from GRIB code. 55C INTFAU - Prepare to interpolate unpacked input field. 56C INTFBU - Interpolate unpacked input field. 57C INTLOG - Log error message. 58C MKFRAME - Create a 'frame' from a rectangular field. 59C 60C 61C Author 62C ------ 63C 64C J.D.Chambers ECMWF February 2001 65C 66C 67C----< 68C -----------------------------------------------------------------| 69C 70 IMPLICIT NONE 71C 72#include "parim.h" 73#include "nifld.common" 74#include "nofld.common" 75#include "intf.h" 76#include "current.h" 77C 78C Parameters 79C 80 INTEGER JPROUTINE 81 PARAMETER (JPROUTINE = 40170 ) 82C 83C Function arguments 84C 85 INTEGER KNVALS, KUGRIB(*), KVGRIB(*), OUTLENU, OUTLENV 86 REAL PUVSH(KNVALS*2), ZNFLDO(*) 87C 88C Local variables 89C 90 CHARACTER*1 HOLDTYP 91 CHARACTER*1 HTYPE 92 INTEGER ILENF 93 INTEGER IRET 94 INTEGER ISIZE 95 INTEGER NCOUNT 96 INTEGER NGAUSS 97 INTEGER NLAT 98 INTEGER NLON 99 INTEGER NTRUNC 100 INTEGER NUMPTS 101 INTEGER NUVFLAG 102 LOGICAL LFRAME 103 LOGICAL LOLDWIND 104 REAL AREA(4) 105 REAL EAST 106 REAL EW 107 REAL GRID(2) 108 REAL NS 109 REAL POLE(2) 110 REAL WEST 111C 112 REAL RGGRID, SWORK 113 POINTER (IRGGRID, RGGRID(1) ) 114 POINTER (ISWORK, SWORK(1) ) 115C 116 INTEGER KPTS(JPGTRUNC*2) 117 REAL GLATS(JPGTRUNC*2) 118C 119 DATA IRGGRID/0/, ISWORK/0/ 120 SAVE IRGGRID, ISWORK 121C 122C Externals 123C 124 INTEGER INTFAU, INTFBU, INTUVDH, HSH2GG !, FIXAREA 125 INTEGER HIRLAMW 126C 127C -----------------------------------------------------------------| 128C* Section 1. Initialise. 129C -----------------------------------------------------------------| 130C 131 100 CONTINUE 132C 133 INTUVXH = 0 134C 135 LFRAME = LNOFRAME.AND. 136 X ((NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPGAUSSIAN).OR. 137 X (NOREPR.EQ.JPREGROT ).OR.(NOREPR.EQ.JPFGGROT ) ) 138C 139 IF( LNOROTA ) GOTO 300 140C 141C -----------------------------------------------------------------| 142C* Section 2. Spectral to grid-point with no rotation 143C -----------------------------------------------------------------| 144C 145 200 CONTINUE 146C 147 CALL INTLOG(JP_DEBUG, 148 X 'INTUVXH: Interoplate U & V fields with no rotation',JPQUIET) 149C 150C Unpack and interpolate U field 151C 152 NIFORM = 0 153 NIPARAM = JP_U 154 LWIND = .TRUE. 155 LOLDWIND = LWINDSET 156 LWINDSET = .TRUE. 157C 158 IRET = INTFAU( PUVSH, KNVALS) 159 IF( IRET.NE.0 ) THEN 160 CALL INTLOG(JP_ERROR, 161 X 'INTUVXH: Prepare to interpolate failed.',JPQUIET) 162 INTUVXH = JPROUTINE + 2 163 GOTO 900 164 ENDIF 165C 166 IRET = INTFBU( PUVSH, KNVALS, ZNFLDO, ILENF) 167 IF( IRET.NE.0 ) THEN 168 CALL INTLOG(JP_ERROR,'INTUVXH: Interpolation failed.',JPQUIET) 169 INTUVXH = JPROUTINE + 2 170 GOTO 900 171 ENDIF 172C 173C Unpack and interpolate V field 174C 175 NIPARAM = JP_V 176C 177 IRET = INTFAU( PUVSH(1+KNVALS), KNVALS) 178 IF( IRET.NE.0 ) THEN 179 CALL INTLOG(JP_ERROR, 180 X 'INTUVXH: Prepare to interpolate failed.',JPQUIET) 181 INTUVXH = JPROUTINE + 2 182 GOTO 900 183 ENDIF 184C 185 IRET = INTFBU( PUVSH(1+KNVALS), KNVALS, ZNFLDO(1+ILENF), ILENF) 186 IF( IRET.NE.0 ) THEN 187 CALL INTLOG(JP_ERROR,'INTUVXH: Interpolation failed.',JPQUIET) 188 INTUVXH = JPROUTINE + 2 189 GOTO 900 190 ENDIF 191C 192C Reset the input format flag 193C 194 NIFORM = 1 195C 196C Get some scratch memory for the U fields before interpolation 197C (to prevent overwriting V in ZNFLDO during interpolation) 198C 199cs ISIZE = ILENF 200cs CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET) 201cs IF( IRET.NE.0 ) THEN 202cs CALL INTLOG(JP_ERROR, 203cs X 'INTUVXH: memory alloc for scratch memory failed',JPQUIET) 204cs INTUVXH = JPROUTINE + 2 205cs GOTO 900 206cs ENDIF 207C 208 LWINDSET = LOLDWIND 209C 210C 211C If a 'frame' has been specified, build the frame 212C 213 IF( LFRAME ) THEN 214 LIMISSV = .TRUE. 215 NLON = 1 + NINT(FLOAT(NOAREA(JPEAST) - NOAREA(JPWEST)) / 216 X NOGRID(JPWESTEP)) 217 NLAT = 1 + NINT(FLOAT(NOAREA(JPNORTH) - NOAREA(JPSOUTH)) / 218 X NOGRID(JPNSSTEP)) 219 CALL MKFRAME(NLON,NLAT,ZNFLDO,RMISSGV,NOFRAME) 220 CALL MKFRAME(NLON,NLAT,ZNFLDO(1+ILENF),RMISSGV,NOFRAME) 221 ENDIF 222 223C Code data into GRIB 224cs DO LOOP = 1, ILENF 225cs SWORK(LOOP) = ZNFLDO(LOOP) 226cs ENDDO 227C 228cs IRET = INTUVDH(SWORK,ILENF,KUGRIB,OUTLENU,'C',JP_U) 229 IRET = INTUVDH(ZNFLDO,ILENF,KUGRIB,OUTLENU,'C',JP_U) 230 IF( IRET.NE.0 ) THEN 231 CALL INTLOG(JP_ERROR, 232 X 'INTUVXH: Wind component into GRIB encoding fail',IRET) 233 INTUVXH = JPROUTINE + 2 234 GOTO 900 235 ENDIF 236C 237 IRET = INTUVDH(ZNFLDO(1+ILENF),ILENF,KVGRIB,OUTLENV,'C',JP_V) 238 IF( IRET.NE.0 ) THEN 239 CALL INTLOG(JP_ERROR, 240 X 'INTUVXH: Wind component into GRIB encoding fail',IRET) 241 INTUVXH = JPROUTINE + 2 242 GOTO 900 243 ENDIF 244C 245 GOTO 900 246C 247C -----------------------------------------------------------------| 248C* Section 3. Initialise spectral to grid-point with rotation 249C -----------------------------------------------------------------| 250C 251 300 CONTINUE 252C 253 IF( .NOT.LUSEHIR ) THEN 254 CALL INTLOG(JP_ERROR, 255 X 'INTUVXH : Unable to rotate spectral U or V:',JPQUIET) 256 INTUVXH = JPROUTINE + 3 257 GOTO 900 258 ENDIF 259C 260 IF( (NOREPR.NE.JPREGROT).AND.(NOREPR.NE.JPREGULAR) ) THEN 261 CALL INTLOG(JP_ERROR, 262 X 'INTUVXH : For U/V, only regular lat/long',JPQUIET) 263 CALL INTLOG(JP_ERROR, 264 X 'INTUVXH : output rotated grids allowed',JPQUIET) 265 INTUVXH = JPROUTINE + 3 266 GOTO 900 267 ENDIF 268C 269 CALL INTLOG(JP_DEBUG,'INTUVXH: Rotate the U & V fields',JPQUIET) 270 CALL INTLOG(JP_DEBUG,'INTUVXH: South pole lat ',NOROTA(1)) 271 CALL INTLOG(JP_DEBUG,'INTUVXH: South pole long ',NOROTA(2)) 272C 273C Fill area limits (handles case when default 0/0/0/0 given) 274C 275cs IRET = FIXAREA() 276cs IF( IRET.NE.0 ) THEN 277cs CALL INTLOG(JP_ERROR,'INTUVXH: area fixup failed',JPQUIET) 278cs INTUVXH = JPROUTINE + 3 279cs GOTO 900 280cs ENDIF 281C 282 AREA(1) = REAL(NOAREA(1))/PPMULT 283 AREA(2) = REAL(NOAREA(2))/PPMULT 284 AREA(3) = REAL(NOAREA(3))/PPMULT 285 AREA(4) = REAL(NOAREA(4))/PPMULT 286C 287 GRID(1) = REAL(NOGRID(1))/PPMULT 288 GRID(2) = REAL(NOGRID(2))/PPMULT 289C 290 POLE(1) = REAL(NOROTA(1))/PPMULT 291 POLE(2) = REAL(NOROTA(2))/PPMULT 292C 293C -----------------------------------------------------------------| 294C* Section 4. Convert spectral to suitable global reduced gaussian 295C -----------------------------------------------------------------| 296C 297 400 CONTINUE 298C 299 NTRUNC = NIRESO 300 NGAUSS = 0 301 HTYPE = '' 302 NS = 0. 303 EW = 0. 304 IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,GLATS,ISIZE) 305 IF( IRET.NE.0 ) THEN 306 CALL INTLOG(JP_ERROR, 307 X 'INTUVXH: problem getting data for reduced grid',NTRUNC) 308 INTUVXH = JPROUTINE + 4 309 GOTO 900 310 ENDIF 311 NCOUNT = ISIZE 312C 313C Dynamically allocate memory for global reduced gaussian grid 314C 315 CALL JMEMHAN( 18, IRGGRID, (NCOUNT*2), 1, IRET) 316 IF( IRET.NE.0 ) THEN 317 CALL INTLOG(JP_ERROR, 318 X 'INTUVXH: memory alloc for reduced grid fail',JPQUIET) 319 INTUVXH = JPROUTINE + 4 320 GOTO 900 321 ENDIF 322C 323C Set flag to show field is a wind component 324C 325 NUVFLAG = 1 326C 327C Create the reduced gaussian grid 328C 329 HOLDTYP = HOGAUST 330 WEST = 0.0 331 EAST = 360.0 - (360.0/FLOAT(KPTS(NGAUSS))) 332C 333C U component 334C 335 CALL JAGGGP(PUVSH,NTRUNC,GLATS(1),GLATS(NGAUSS*2),WEST, 336 X EAST,NGAUSS,HTYPE,KPTS,RGGRID,NUVFLAG,IRET) 337 IF( IRET.NE.0 ) THEN 338 CALL INTLOG(JP_ERROR, 339 X 'INTUVXH: spectral to reduced gaussian failed',JPQUIET) 340 INTUVXH = JPROUTINE + 4 341 GOTO 900 342 ENDIF 343C 344 HOGAUST = HOLDTYP 345C 346C V component 347C 348 CALL JAGGGP(PUVSH(1+KNVALS),NTRUNC,GLATS(1),GLATS(NGAUSS*2),WEST, 349 X EAST,NGAUSS,HTYPE,KPTS,RGGRID(1+NCOUNT),NUVFLAG,IRET) 350 IF( IRET.NE.0 ) THEN 351 CALL INTLOG(JP_ERROR, 352 X 'INTUVXH: spectral to reduced gaussian failed',JPQUIET) 353 INTUVXH = JPROUTINE + 4 354 GOTO 900 355 ENDIF 356C 357 HOGAUST = HOLDTYP 358C 359C -----------------------------------------------------------------| 360C* Section 5. Rotate using 12-point horizontal interpolation 361C -----------------------------------------------------------------| 362C 363 500 CONTINUE 364C 365C Dynamically allocate memory for rotated lat/long grid 366C 367 NLON = 1 + NINT(FLOAT(NOAREA(JPEAST) - NOAREA(JPWEST)) / 368 X NOGRID(JPWESTEP)) 369 NLAT = 1 + NINT(FLOAT(NOAREA(JPNORTH) - NOAREA(JPSOUTH)) / 370 X NOGRID(JPNSSTEP)) 371C 372 NUMPTS = NLON * NLAT 373 ISIZE = NUMPTS * 2 374 CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET) 375 IF( IRET.NE.0 ) THEN 376 CALL INTLOG(JP_ERROR, 377 X 'INTUVXH: memory alloc for lat/long grid fail',JPQUIET) 378 INTUVXH = JPROUTINE + 5 379 GOTO 900 380 ENDIF 381C 382 IRET = HIRLAMW(LO12PT,RGGRID,RGGRID(1+NCOUNT),NCOUNT,NGAUSS,HTYPE, 383 X AREA,POLE,GRID,SWORK,SWORK(1+NUMPTS),NUMPTS,NLON, 384 X NLAT) 385C 386 IF( IRET.NE.0 ) THEN 387 CALL INTLOG(JP_ERROR, 388 X 'INTUVXH: HIRLAMW rotation failed',JPQUIET) 389 INTUVXH = JPROUTINE + 5 390 GOTO 900 391 ENDIF 392C 393C -----------------------------------------------------------------| 394C* Section 6. Pack the fields into GRIB format 395C -----------------------------------------------------------------| 396C 397 600 CONTINUE 398C 399C Reset the input format flag 400C 401 NIFORM = 1 402C 403C Set the components flag for rotated U and V coefficients 404C 405 ISEC2(19) = 8 406C 407C If a 'frame' has been specified, build the frame 408C 409 IF( LFRAME ) THEN 410 LIMISSV = .TRUE. 411 CALL MKFRAME(NLON,NLAT,SWORK,RMISSGV,NOFRAME) 412 CALL MKFRAME(NLON,NLAT,SWORK(1+NUMPTS),RMISSGV,NOFRAME) 413 ENDIF 414C 415 IRET = INTUVDH(SWORK,NUMPTS,KUGRIB,OUTLENU,'C',JP_U) 416 IF( IRET.NE.0 ) THEN 417 CALL INTLOG(JP_ERROR, 418 X 'INTUVXH: Wind component into GRIB encoding fail',IRET) 419 INTUVXH = JPROUTINE + 6 420 GOTO 900 421 ENDIF 422C 423 IRET = INTUVDH(SWORK(1+NUMPTS),NUMPTS,KVGRIB,OUTLENV,'C',JP_V) 424 IF( IRET.NE.0 ) THEN 425 CALL INTLOG(JP_ERROR, 426 X 'INTUVXH: Wind component into GRIB encoding fail',IRET) 427 INTUVXH = JPROUTINE + 6 428 GOTO 900 429 ENDIF 430C 431C -----------------------------------------------------------------| 432C* Section 9. Return 433C -----------------------------------------------------------------| 434C 435 900 CONTINUE 436C 437 RETURN 438 END 439