1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
INTEGER FUNCTION INTUVGH(KUGRIBI, KVGRIBI, INLEN,
X KUGRIBO, KVGRIBO,
X OUTLENU, OUTLENV)
C
C---->
C**** INTUVGH
C
C Purpose
C -------
C
C Interpolate GRIB format input U and V component reduced gaussian
C fields to rotated GRIB format U and V fields.
C
C
C Interface
C ---------
C
C IRET = INTUVGH(KUGRIBI,KVGRIBI,INLEN,KUGRIBO,KVGRIBO,OUTLENU,OUTLENV)
C
C Input
C -----
C
C KUGRIBI - Input U component field (reduced gaussian, GRIB format).
C KVGRIBI - Input V component field (reduced gaussian, GRIB format).
C INLEN - Input field length (words).
C
C
C Output
C ------
C
C KUGRIBO - Output U field (GRIB format).
C KVGRIBO - Output V field (GRIB format).
C OUTLENU - Output U field length (words).
C OUTLENV - Output V field length (words).
C
C
C Method
C ------
C
C Use HIRLAM rotation scheme.
C
C Note that a common block is used in intf.h to hold the U/V
C fields before interpolation.
C
C Externals
C ---------
C
C IBASINI - Ensure basic interpolation setup is done.
C INTUVDH - Encodes/decodes data into/from GRIB code.
C JMEMHAN - Allocate scratch memory.
C GRIBEX - GRIB decoding/encoding.
C ISCRSZ - Calculate number of values in generated field.
C FIXAREA - Fixup area definition to correspond to grid definitions
C INTLOG - Log error message.
C MKFRAME - Create a 'frame' from a rectangular field.
C RESET_C - Reset interpolation handling options using GRIB product.
C INSANE - Ensure no outrageous values given for interpolation.
C HRG2GGW - Creates rotated gaussian U and V from reduced gaussian
C HLL2LLW - Creates rotated lat/long U and V from lat/long
C
C
C Author
C ------
C
C J.D.Chambers ECMWF October 2001
C
C
C----<
C -----------------------------------------------------------------|
C
IMPLICIT NONE
C
C Function arguments
C
INTEGER KUGRIBI(*), KVGRIBI(*), INLEN
INTEGER KUGRIBO(*), KVGRIBO(*), OUTLENU, OUTLENV
C
#include "parim.h"
#include "nifld.common"
#include "nofld.common"
#include "grfixed.h"
#include "intf.h"
#include "intlog.h"
#include "current.h"
C
C Parameters
C
INTEGER JPROUTINE, JPALLOC, JPDEALL, JPSCR3, JPSCR4,JPSCR5
PARAMETER (JPROUTINE = 40170 )
PARAMETER (JPALLOC = 1)
PARAMETER (JPDEALL = 0)
PARAMETER (JPSCR3 = 3)
PARAMETER (JPSCR4 = 4)
PARAMETER (JPSCR5 = 5)
C
C Local variables
C
CHARACTER*1 HFUNC
REAL OLDGRID(2), AREA(4), GRID(2), POLE(2)
LOGICAL LOLDWIND, LFRAME
INTEGER IRET,KPR,ISZOUV,ISIZE,ISAME,IPARAM
INTEGER LOOP, NGAUSS
INTEGER IPINU, IPINV, NUMPTS, NLON, NLAT, IN_U, IN_V
#ifdef POINTER_64
INTEGER*8 IOLDUV, ISWORK
#endif
REAL OLDUV, SWORK
POINTER ( IOUV, OLDUV(1) )
POINTER (ISWORK, SWORK(1) )
C
C Externals
C
INTEGER RESET_C,ISCRSZ,FIXAREA,IBASINI,INSANE,INTUVDH,HIRLAMW
INTEGER HRG2GGW, HLL2LLW
C
C -----------------------------------------------------------------|
C* Section 1. Initialise
C -----------------------------------------------------------------|
C
100 CONTINUE
INTUVGH = 0
IRET = 0
KPR = 0
C
LFRAME = LNOFRAME.AND.
X ((NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPREGROT))
C
C Ensure that basic initialisation has been done
C
IRET = IBASINI(0)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'INTUVGH: basic initialise failed',JPQUIET)
INTUVGH = IRET
GOTO 900
ENDIF
C
C -----------------------------------------------------------------|
C* Section 2. Unpack the U/V fields.
C -----------------------------------------------------------------|
C
200 CONTINUE
C
C Need to establish input characteristics, so unpack GRIB sections 1
C and 2.
C
IPARAM = 0
IPINU = 1
IRET = INTUVDH(OLDUV(IPINU),ISZOUV,KUGRIBI,INLEN,'I',IPARAM)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'INTUVGH: GRIB header decode failed',IRET)
INTUVGH = IRET
GOTO 900
ENDIF
C
C Check whether the input field is an ECMWF reduced gaussian or
C lat/long U field
C
IF( (ISEC1(1).NE.128).OR.
X ((ISEC1(6).NE.JP_U).AND.(ISEC1(6).NE.JP_10U)) ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Input is not ECMWF U field.',JPQUIET)
INTUVGH = JPROUTINE + 2
GOTO 900
ENDIF
IF( (ISEC2(1).NE.4).AND.
X ((ISEC2(17).NE.0) .OR.(ISEC2(1).NE.0)) ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Input not gaussian or lat/long.',JPQUIET)
INTUVGH = JPROUTINE + 2
GOTO 900
ENDIF
C
C Get scratch memory for U/V unpacked fields.
C U/V memory areas are adjacent.
C
IF( ISEC2(17).EQ.1 ) THEN
ISZOUV = 0
DO LOOP = 1, ISEC2(3)
ISZOUV = ISZOUV + ISEC2(22+LOOP)
ENDDO
ELSE
ISZOUV = ISEC2(2) * ISEC2(3)
ENDIF
C
IPINU = 1
IPINV = 1 + ISZOUV
CALL JMEMHAN( JPSCR4, IOUV, ISZOUV*2, JPALLOC, IRET)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Scratch memory type 4 allocation failed.',JPQUIET)
INTUVGH = IRET
GOTO 900
ENDIF
C
C Decode input U and V
C
IRET = INTUVDH(OLDUV(IPINU),ISZOUV,KUGRIBI,INLEN,'D',IPARAM)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'INTUVGH: Vorticity decoding failed',IRET)
INTUVGH = IRET
GOTO 900
ENDIF
IN_U = ISEC1(6)
C
IRET = INTUVDH(OLDUV(IPINV),ISZOUV,KVGRIBI,INLEN,'D',IPARAM)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'INTUVGH: Divergence decoding failed',IRET)
INTUVGH = IRET
GOTO 900
ENDIF
IN_V = ISEC1(6)
C
C Setup interpolation options from input GRIB characteristics.
C
IRET = RESET_C( ISEC1, ISEC2, ZSEC2, ISEC4)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Setup interp. options from GRIB failed.',JPQUIET)
INTUVGH = IRET
GOTO 900
ENDIF
C
C Check that no outrageous values given for interpolation
C
ISAME = INSANE()
IF( ISAME.GT.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Interpolation cannot use given values.',JPQUIET)
INTUVGH = ISAME
GOTO 900
ENDIF
C
IRET = FIXAREA()
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'INTUVGH: Fixarea failed.',JPQUIET)
INTUVGH = IRET
GOTO 900
ENDIF
C
LWIND = .TRUE.
LOLDWIND = LWINDSET
LWINDSET = .TRUE.
C
DO LOOP = 1, 4
AREA(LOOP) = REAL(NOAREA(LOOP))/PPMULT
ENDDO
DO LOOP = 1, 2
GRID(LOOP) = REAL(NOGRID(LOOP))/PPMULT
POLE(LOOP) = REAL(NOROTA(LOOP))/PPMULT
ENDDO
C
C -----------------------------------------------------------------|
C* Section 3. Generate interpolated lat/long U and V fields.
C -----------------------------------------------------------------|
C
300 CONTINUE
C
IF( (NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPREGROT) ) THEN
C
C Dynamically allocate scrath space for rotated lat/long grid
C
NLON = 1 + NINT(FLOAT(NOAREA(JPEAST) - NOAREA(JPWEST)) /
X NOGRID(JPWESTEP))
NLAT = 1 + NINT(FLOAT(NOAREA(JPNORTH) - NOAREA(JPSOUTH)) /
X NOGRID(JPNSSTEP))
C
NOWE = NLON
NONS = NLAT
NUMPTS = NLON * NLAT
ISIZE = NUMPTS * 2
CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: memory alloc for lat/long grid fail',JPQUIET)
INTUVGH = JPROUTINE + 3
GOTO 900
ENDIF
C
C Rotate reduced gaussian to lat/long
C
IF( (NIREPR.EQ.JPGAUSSIAN).OR.(NIREPR.EQ.JPQUASI) ) THEN
IRET = HIRLAMW(LO12PT,
X OLDUV(IPINU),OLDUV(IPINV),ISZOUV,
X NOGAUSS,AREA,POLE,GRID,
X SWORK,SWORK(1+NUMPTS),NUMPTS,NLON,NLAT)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: HIRLAMW rotation failed',JPQUIET)
INTUVGH = JPROUTINE + 3
GOTO 900
ENDIF
C
C Rotate lat/long to lat/long
C
ELSE
OLDGRID(1) = REAL(NIGRID(1)) / PPMULT
OLDGRID(2) = REAL(NIGRID(2)) / PPMULT
IRET = HLL2LLW(LO12PT,OLDUV(IPINU),OLDUV(IPINV),
X OLDGRID,AREA,POLE,GRID,
X SWORK,SWORK(1+NUMPTS),NUMPTS,NLON,NLAT)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: HLL2LLW rotation failed',JPQUIET)
INTUVGH = JPROUTINE + 3
GOTO 900
ENDIF
ENDIF
C
C If a 'frame' has been specified, build the frame
C
IF( LFRAME ) THEN
LIMISSV = .TRUE.
CALL MKFRAME(NLON,NLAT,SWORK,RMISSGV,NOFRAME)
CALL MKFRAME(NLON,NLAT,SWORK(1+NUMPTS),RMISSGV,NOFRAME)
ENDIF
C
ELSE
C
C -----------------------------------------------------------------|
C* Section 4. Generate interpolated gaussian U and V fields.
C -----------------------------------------------------------------|
C
400 CONTINUE
C
IF( LDEBUG ) CALL INTLOG(JP_DEBUG,
X 'INTUVGH: Gaussian to gaussian',JPQUIET)
C
C Dynamically allocate memory for rotated gaussian grids
C
NUMPTS = NOGAUSS * NOGAUSS
ISIZE = 2 * NUMPTS * 8
CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET)
IF( IRET.NE.0 ) THEN
IF( LDEBUG ) CALL INTLOG(JP_ERROR,
X 'INTUVGH: memory alloc for gaussian grid fail',JPQUIET)
INTUVGH = JPROUTINE + 4
GOTO 900
ENDIF
C
NGAUSS = ISEC2(10)
IRET = HRG2GGW(LO12PT,
X OLDUV(IPINU),OLDUV(IPINV),ISZOUV,
X NGAUSS,AREA,POLE,NOGAUSS,HOGAUST,
X SWORK,SWORK(1+ISIZE/2),ISIZE,NUMPTS)
IF( IRET.NE.0 ) THEN
IF( LDEBUG ) CALL INTLOG(JP_ERROR,
X 'INTUVGH: HRG2GGW rotation failed',JPQUIET)
INTUVGH = JPROUTINE + 4
GOTO 900
ENDIF
C
ENDIF
C
C -----------------------------------------------------------------|
C* Section 5. Pack the fields into GRIB format
C -----------------------------------------------------------------|
C
500 CONTINUE
C
IF( ((NOREPR.EQ.JPQUASI).OR.
X (NOREPR.EQ.JPQGGROT).OR.
X (NOREPR.EQ.JPFGGROT).OR.
X (NOREPR.EQ.JPGAUSSIAN)).AND.
X ( (NOROTA(1).NE.9000000).AND.(NOROTA(2).NE.0) ) ) THEN
ISEC2(1) = JPFGGROT
ELSE
ISEC2(1) = NOREPR
ENDIF
C
ISEC4(1) = NUMPTS
C
C Reset the input format flag
C
NIFORM = 1
C
C Set the components flag for rotated U and V coefficients
C
ISEC2(19) = 8
C
IRET = INTUVDH(SWORK,NUMPTS,KUGRIBO,OUTLENU,'C',IN_U)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Wind component into GRIB encoding fail',IRET)
INTUVGH = JPROUTINE + 4
GOTO 900
ENDIF
C
IRET = INTUVDH(SWORK(1+(ISIZE/2)),NUMPTS,KVGRIBO,OUTLENV,'C',IN_V)
IF( IRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,
X 'INTUVGH: Wind component into GRIB encoding fail',IRET)
INTUVGH = JPROUTINE + 4
GOTO 900
ENDIF
C
C -----------------------------------------------------------------|
C* Section 9. Return
C -----------------------------------------------------------------|
C
900 CONTINUE
C
C Clear change flags for next product processing and reset wind flag
C
LCHANGE = .FALSE.
LSMCHNG = .FALSE.
LWINDSET = LOLDWIND
C
RETURN
END
|
