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