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 INTWAVE2(INGRIB,INLEN,OUTGRIB,OUTLEN)
12C
13C---->
14C**** INTWAVE2
15C
16C     Purpose
17C     -------
18C
19C     Interpolate an ECMWF wave field wave field or a reduced
20C     latitude-longitude field
21C
22C
23C     Interface
24C     ---------
25C
26C     IRET = INTWAVE2(INGRIB,INLEN,OUTGRIB,OUTLEN)
27C
28C     Input
29C     -----
30C
31C     INGRIB - Input field (packed).
32C     INLEN  - Input field length (words).
33C
34C
35C     Output
36C     ------
37C
38C     OUTGRIB - Output field (packed).
39C     OUTLEN  - Output field length (words).
40C
41C
42C     Method
43C     ------
44C
45C     Call interpolation routine; then repack into GRIB format.
46C
47C
48C     Externals
49C     ---------
50C
51C     FIXAREA - Fixup area definition to correspond to grid definitions
52C     GRIBEX  - Decode/encode GRIB product.
53C     JMEMHAN - Handles memory allocation.
54C     JDEBUG  - Checks environment to switch on/off debug
55C     INTLOG  - Log error message.
56C     WAVEXX1 - Interpolate wave fields (except 2D spectra)
57C     WV2DXXX - Interpolate wave 2D spectra fields
58C     MKFRAME - Create a 'frame' from a rectangular field
59C     MKBITMP - Apply a bitmap to a rectangular field
60C
61C
62C     Author
63C     ------
64C
65C     S.Curic     ECMWF     Jun 2009
66C
67C
68C----<
69C
70      IMPLICIT NONE
71C
72C     Function arguments
73      INTEGER INLEN,OUTLEN
74      INTEGER INGRIB(*),OUTGRIB(*)
75C
76#include "parim.h"
77#include "nifld.common"
78#include "nofld.common"
79#include "intf.h"
80#include "grfixed.h"
81#include "jparams.h"
82C
83C     Parameters
84      INTEGER JPROUTINE, JP2DSP, JP2DSPQ, JPMAXLT
85
86      PARAMETER (JPROUTINE = 40200 )
87      PARAMETER (JP2DSP  = 250 )  ! for 2D wave spectra (whole)
88      PARAMETER (JP2DSPQ = 251 )  ! for 2D wave spectra (single)
89      PARAMETER (JPMAXLT = 1801 )
90C
91C     Local variables
92C
93      CHARACTER*1 HFUNC
94
95      REAL ZMISS
96      DATA ZMISS/-9999999.0/
97      INTEGER IWORD, IRET, ISIZE, NSPEC, LAT1, LOOP, IERR
98      INTEGER KNUM, NUM_E_W, NUM_N_S, NFULLNS, NFULLEW
99      REAL GRIDWE, GRIDNS, NORTH, SOUTH, WEST, EAST
100      INTEGER NUMPTS
101      DIMENSION NUMPTS(JPMAXLT)
102#ifndef _CRAYFTN
103#ifdef POINTER_64
104      INTEGER*8 INEWAVE, IDISTNW, IZ2DSP, INEWIDX
105#endif
106#endif
107      REAL NEWAVE
108      POINTER ( INEWAVE, NEWAVE )
109      DIMENSION NEWAVE( 1 )
110      REAL*4 DISTNEW
111      POINTER ( IDISTNW, DISTNEW )
112      DIMENSION DISTNEW( 1 )
113      REAL Z2DSP
114      POINTER ( IZ2DSP, Z2DSP )
115      DIMENSION Z2DSP( 1 )
116      INTEGER NEWIDX
117      POINTER ( INEWIDX, NEWIDX )
118      DIMENSION NEWIDX( 1 )
119      REAL RNLAT, RSLAT
120      INTEGER NUMTABL, NUMPROD
121      LOGICAL LCOEFFS, L98WAVE, LFRAME, LBITMP
122      LOGICAL LDEBUG
123      INTEGER NPARAM
124      INTEGER MISSLAT, KOLDNUM
125C
126C     Externals
127      INTEGER WAVEXX2, WV2DXX2, AREACHK, FIXAREA
128      INTEGER NUMPTNS, NUMPTWE, JNORSGG, MKBITMP
129C
130C ------------------------------------------------------------------
131C     Section 1.   Initialise
132C ------------------------------------------------------------------
133C
134  100 CONTINUE
135C
136      INTWAVE2 = 0
137C
138      LCOEFFS = .FALSE.
139C
140C     Check if debug option turned on
141      CALL JDEBUG()
142      LDEBUG = NDBG.GT.0
143C
144C     Allocate work array ZNFELDI if not already done.
145C
146      IF( IZNJDCI.NE.1952999238 ) THEN
147        CALL JMEMHAN( 19, IZNFLDI, JPEXPAND, 1, IRET)
148        IF( IRET.NE.0 ) THEN
149          CALL INTLOG(JP_WARN,'INTWAVE2: ZNFELDI allocate fail',JPQUIET)
150          INTWAVE2= IRET
151          GOTO 900
152        ENDIF
153        IZNJDCI = 1952999238
154      ENDIF
155C
156C     Unpack GRIB message headers (using ZMISS as missing data value).
157C
158cs      IWORD   = INLEN
159      IRET = 1
160      ISEC3(2) = NINT(ZMISS)
161      ZSEC3(2) = ZMISS
162      CALL GRIBEX (ISEC0,ISEC1,ISEC2,ZSEC2,ISEC3,ZSEC3,ISEC4,
163     X             ZNFELDI,JPEXPAND,INGRIB,INLEN,IWORD,'D',IRET)
164
165      IF(IRET.GT.0) THEN
166       CALL INTLOG(JP_ERROR,'INTWAVE2: GRIBEX decoding failed.',JPQUIET)
167        INTWAVE2 = JPROUTINE + 1
168        GOTO 900
169      ENDIF
170
171      NPARAM = ISEC1(6)
172C
173C     Check that the field can be handled
174C
175      NUMTABL = ISEC1(2)*1000 + ISEC1(1)
176      NUMPROD = NUMTABL*1000 + ISEC1(6)
177      L98WAVE = (NUMTABL.EQ.98140).OR.
178     X          (NUMPROD.EQ.98131229).OR.
179     X          (NUMPROD.EQ.98131232).OR.
180     X          (NIREPR.EQ.26)
181      IF( .NOT.L98WAVE ) THEN
182        CALL INTLOG(JP_WARN,'INTWAVE2: Not an ECMWF wave field',JPQUIET)
183        INTWAVE2 = JPROUTINE + 2
184        GOTO 900
185      ENDIF
186      IF( ISEC2(1).NE.0 ) THEN
187        CALL INTLOG(JP_WARN,
188     X    'INTWAVE2: Not a lat/long field',JPQUIET)
189        INTWAVE2 = JPROUTINE + 3
190        GOTO 900
191      ENDIF
192      IF( (ISEC2(6).NE.128) .OR. (ISEC2(10).LT.10) ) THEN
193        CALL INTLOG(JP_ERROR,
194     X    'INTWAVE2: Cannot handle longitude increment',ISEC2(10))
195        INTWAVE2 = JPROUTINE + 4
196        GOTO 900
197      ENDIF
198C
199C ------------------------------------------------------------------
200C     Section 2.   Sort out interpolation area.
201C ------------------------------------------------------------------
202C
203  200 CONTINUE
204C
205      NSPEC = 1
206      KNUM  = ISEC2(3)
207C
208C     Adjust (sub-)area limits to suit the grid
209C
210      IRET = FIXAREA()
211      IF ( IRET .NE. 0 ) THEN
212        CALL INTLOG(JP_ERROR, 'INTWAVE2: FIXAREA failed.',JPQUIET)
213        INTWAVE2 = JPROUTINE + 5
214        GOTO 900
215      ENDIF
216      NORTH = FLOAT(NOAREA(1))/PPMULT
217      WEST  = FLOAT(NOAREA(2))/PPMULT
218      SOUTH = FLOAT(NOAREA(3))/PPMULT
219      EAST  = FLOAT(NOAREA(4))/PPMULT
220C
221C     Calculate the number of points E-W and N-S in new grid area
222C
223      IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN)) THEN
224        GRIDWE  = NOGAUSS
225        GRIDNS  = 0
226      ELSE IF( NOREPR.EQ.JPREDLL ) THEN
227        GRIDWE  = 360.0/NOLPTS((NOREDLL/2)+1)
228        GRIDNS  = FLOAT(NOGRID(2))/PPMULT
229      ELSE
230        GRIDWE  = FLOAT(NOGRID(1))/PPMULT
231        GRIDNS  = FLOAT(NOGRID(2))/PPMULT
232      ENDIF
233C
234      IRET = AREACHK(GRIDWE,GRIDNS,NORTH,WEST,SOUTH,EAST)
235      IF( IRET.NE.0 ) THEN
236        INTWAVE2 = JPROUTINE + 6
237        GOTO 900
238      ENDIF
239      NOAREA(1) = NINT(NORTH*PPMULT)
240      NOAREA(2) = NINT(WEST*PPMULT)
241      NOAREA(3) = NINT(SOUTH*PPMULT)
242      NOAREA(4) = NINT(EAST*PPMULT)
243C
244      IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN)) THEN
245        GRIDWE  = 360.0/NOLPTS(NOGAUSS)
246        NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE)
247        NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1
248        NUM_N_S = JNORSGG(SOUTH,ROGAUSS,NOGAUSS,1) -
249     X            JNORSGG(NORTH,ROGAUSS,NOGAUSS,0) + 1
250        NFULLNS = NOGAUSS*2
251      ELSE IF( NOREPR.EQ.JPREDLL ) THEN
252        NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE)
253        NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1
254        NUM_N_S = NUMPTNS(NORTH,SOUTH,GRIDNS)
255        NFULLNS = NOREDLL
256      ELSE
257        NUM_E_W = NUMPTWE(WEST,EAST,GRIDWE)
258        NFULLEW = NUMPTWE(0.0,360.0,GRIDWE) - 1
259        NUM_N_S = NUMPTNS(NORTH,SOUTH,GRIDNS)
260        NFULLNS = NUM_N_S
261      ENDIF
262C
263C
264C--------------------------------------------------------------------
265        KNUM = NINS
266cs        KNUM = ISEC2(3)
267      IF( ISEC2(17).EQ.1 ) THEN
268C
269C       Input field is a reduced latitude/longitude grid
270C
271C       .. but it may be 'pseudo-gaussian' in layout
272C       (ie global, symmetric about the equator but no latitude
273C        at the equator)
274C
275        IF( (ISEC2(4).NE.90000).AND.(MOD(ISEC2(3),2).EQ.0) ) THEN
276C
277          DO LOOP = 1, ISEC2(3)
278            NUMPTS(LOOP) = ISEC2(22+LOOP)
279          ENDDO
280C
281        ELSE
282C
283          MISSLAT = (90000 - ISEC2(4))/ISEC2(10)
284          DO LOOP = 1, MISSLAT
285            NUMPTS(LOOP)    = 0
286          ENDDO
287          KOLDNUM = 1 + (90000 - ISEC2(7))/ISEC2(10)
288          DO LOOP = 1, (KOLDNUM-MISSLAT)
289            NUMPTS(LOOP+MISSLAT) = ISEC2(22+LOOP)
290          ENDDO
291          DO LOOP = (KOLDNUM+1), KNUM
292            NUMPTS(LOOP)    = 0
293          ENDDO
294        ENDIF
295      ENDIF
296C
297C
298C ------------------------------------------------------------------
299C     Section 3.   Interpolate wave fields other than 2D spectra.
300C ------------------------------------------------------------------
301C
302  300 CONTINUE
303C
304C     Handle if not 2D spectra ..
305C
306      IF( (ISEC1(6).NE.JP2DSP).AND.(ISEC1(6).NE.JP2DSPQ) ) THEN
307       IF( LDEBUG ) THEN
308       CALL INTLOG(JP_DEBUG,
309     X 'INTWAVE2: Interpolate wave flds other than 2D spectra',JPQUIET)
310       ENDIF
311C
312C       Get some scratch space for the interpolated field.
313C
314        ISIZE = NFULLEW * NUM_N_S
315        CALL JMEMHAN( 3, INEWAVE, ISIZE, 1, IRET)
316        IF ( IRET.NE.0 ) THEN
317          CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET)
318          INTWAVE2 = JPROUTINE + 7
319          GOTO 900
320        ENDIF
321C
322        ISIZE = NFULLEW * NFULLNS * 4
323        CALL JMEMHAN( 4, INEWIDX, ISIZE, 1, IRET)
324        IF ( IRET.NE.0 ) THEN
325          CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET)
326          INTWAVE2 = JPROUTINE + 8
327          GOTO 900
328        ENDIF
329C
330        ISIZE = NFULLEW * NFULLNS * 10
331        CALL JMEMHAN( 5, IDISTNW, ISIZE, 1, IRET)
332        IF ( IRET.NE.0 ) THEN
333          CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET)
334          INTWAVE2 = JPROUTINE + 9
335          GOTO 900
336        ENDIF
337
338C       Interpolate the field
339C
340        IRET = WAVEXX2(NPARAM,KNUM,NUMPTS,NUM_N_S,GRIDNS,GRIDWE,ZNFELDI,
341     X               NEWAVE,NORTH,WEST,ZMISS)
342        IF( IRET.NE.0 ) THEN
343         CALL INTLOG(JP_ERROR,'INTWAVE2: Interpolation failed.',JPQUIET)
344          INTWAVE2 = JPROUTINE + 10
345          GOTO 900
346        ENDIF
347C
348C ------------------------------------------------------------------
349C     Section 4.   Interpolate wave 2D spectra field.
350C ------------------------------------------------------------------
351C
352      ELSE
353C
354  400   CONTINUE
355
356      IF( LDEBUG ) THEN
357      CALL INTLOG(JP_DEBUG,
358     X 'INTWAVE2: Interpolate wave 2D spectra field',JPQUIET)
359      ENDIF
360C
361C       Find number of 2D spectra values at each point
362C
363        IF( ISEC4(8).NE.0 ) LCOEFFS = .TRUE.
364        IF( (ISEC4(8).NE.0).AND.(ISEC1(6).NE.JP2DSPQ) ) THEN
365          NSPEC = ISEC4(50) * ISEC4(51)
366          NIMATR = 1
367        ELSE
368          NSPEC = 1
369          NIMATR = 0
370        ENDIF
371C
372C       Get some scratch space for the interpolated field.
373C
374        ISIZE = NSPEC * NFULLEW * NUM_N_S
375        CALL JMEMHAN( 3, INEWAVE, ISIZE, 1, IRET)
376        IF ( IRET.NE.0 ) THEN
377          CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET)
378          INTWAVE2 = JPROUTINE + 11
379          GOTO 900
380        ENDIF
381C
382        ISIZE = NFULLEW * NFULLNS
383        CALL JMEMHAN( 4, INEWIDX, ISIZE, 1, IRET)
384        IF ( IRET.NE.0 ) THEN
385          CALL INTLOG(JP_ERROR,'INTWAVE2: Get work space fail.',JPQUIET)
386          INTWAVE2 = JPROUTINE + 12
387          GOTO 900
388        ENDIF
389C
390C       Interpolate the field
391C
392        IRET = WV2DXX2(NIMATR,KNUM,NUMPTS,NUM_N_S,GRIDNS,GRIDWE,ZNFELDI,
393     X               NEWAVE,NORTH,WEST,ZMISS)
394        IF( IRET.NE.0 ) THEN
395         CALL INTLOG(JP_ERROR,'INTWAVE2: Interpolation failed.',JPQUIET)
396          INTWAVE2 = JPROUTINE + 13
397          GOTO 900
398        ENDIF
399C
400      ENDIF
401C
402C ------------------------------------------------------------------
403C     Section 5.   Pack the interpolated field into GRIB.
404C ------------------------------------------------------------------
405C
406  500   CONTINUE
407
408C     If a 'bitmap' has been specified, build the bitmap
409C
410      LBITMP = LNOBITMP.AND.
411     X         ((NOREPR.EQ.JPREGROT).OR.(NOREPR.EQ.JPREGULAR).OR.
412     X          (NOREPR.EQ.JPGAUSSIAN))
413
414      IF( LBITMP ) THEN
415        CALL INTLOG(JP_DEBUG,'INTWAVE2: MKBITMP is enabled',JPQUIET)
416        IERR = MKBITMP(NUM_E_W,NUM_N_S,NEWAVE,ZMISS)
417        IF( IERR.NE.0 ) THEN
418          CALL INTLOG(JP_ERROR,'INTFB: Problem applying bitmap',JPQUIET)
419          GOTO 900
420        ENDIF
421      ENDIF
422C     If a 'frame' has been specified, build the frame
423C
424       LFRAME = LNOFRAME.AND.
425     X         ((NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPGAUSSIAN).OR.
426     X          (NOREPR.EQ.JPREGROT ).OR.(NOREPR.EQ.JPFGGROT  ) )
427
428         CALL INTLOG(JP_DEBUG,'INTWAVE2: NUM_E_W ',NUM_E_W)
429         CALL INTLOG(JP_DEBUG,'INTWAVE2: NUM_N_S ',NUM_N_S)
430      IF( LFRAME ) THEN
431           CALL INTLOG(JP_DEBUG,'INTWAVE2: MKFRAME is enabled',JPQUIET)
432           CALL INTLOG(JP_DEBUG,'INTWAVE2: NOFRAME ',NOFRAME)
433         CALL MKFRAME(NUM_E_W,NUM_N_S,NEWAVE,ZMISS,NOFRAME)
434      ENDIF
435C
436C     GRIB sections 2 and 3 (bitmap) included
437C
438      ISEC1(5) = 192
439C
440      IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPREDLL) ) THEN
441C
442C       Reduced gaussian or reduced lat/long with (possible) subarea
443C
444        IF( NOREPR.EQ.JPQUASI ) THEN
445          ISEC2(1) = JPGAUSSIAN
446          ISEC2(10)= NOGAUSS
447        ELSE
448          ISEC2(1) = JPREGULAR
449          ISEC2(10)= NINT(GRIDNS*1000.0)
450        ENDIF
451        ISEC2(2) = 255
452        ISEC2(3) = NUM_N_S
453        ISEC2(5) = NINT(WEST*1000.0)
454        ISEC2(6) = 0
455        ISEC2(8) = NINT(EAST*1000.0)
456        ISEC2(17)= 1
457        ISEC4(1) = 0
458C
459        LAT1 = 0
460 510    CONTINUE
461        LAT1 = LAT1 + 1
462        IF( NOREPR.EQ.JPQUASI ) THEN
463          RNLAT = ROGAUSS(LAT1)
464          RSLAT = ROGAUSS(LAT1+NUM_N_S-1)
465        ELSE
466          RNLAT = 90.0 - (LAT1-1)*GRIDNS
467          RSLAT = 90.0 - (LAT1+NUM_N_S-2)*GRIDNS
468        ENDIF
469        IF( RNLAT.GT.NORTH ) GOTO 510
470C
471        ISEC2(4) = NINT(RNLAT*1000.0)
472        ISEC2(7) = NINT(RSLAT*1000.0)
473C
474        DO LOOP = LAT1, (NUM_N_S+LAT1-1)
475          ISEC2(23+LOOP-LAT1) = NOLPTS(LOOP)
476          ISEC4(1) = ISEC4(1) + ISEC2(23+LOOP-LAT1)
477        ENDDO
478C
479      ELSE
480C
481C       Regular gaussian or lat/long with (possible) subarea
482C
483        IF( NOREPR.EQ.JPGAUSSIAN ) THEN
484          ISEC2(1) = JPGAUSSIAN
485        ELSE
486          ISEC2(1) = JPREGULAR
487        ENDIF
488        ISEC2(2) = NUM_E_W
489        ISEC2(3) = NUM_N_S
490        ISEC2(4) = NINT(NORTH*1000.0)
491        ISEC2(5) = NINT(WEST*1000.0)
492        ISEC2(6) = 128
493        ISEC2(7) = NINT(SOUTH*1000.0)
494        ISEC2(8) = NINT(EAST*1000.0)
495        ISEC2(9) = NINT(GRIDWE*1000.0)
496        IF( NOREPR.EQ.JPGAUSSIAN ) THEN
497          ISEC2(10) = NOGAUSS
498        ELSE
499          ISEC2(10) = NINT(GRIDNS*1000.0)
500        ENDIF
501        ISEC2(17)= 0
502        ISEC4(1) = NUM_E_W * NUM_N_S
503      ENDIF
504C
505C
506      ISEC3(2) = NINT(ZMISS)
507      ZSEC3(2) = ZMISS
508C
509      ISEC4(1) = ISEC4(1) * NSPEC
510      ISEC4(2) = NOACC
511      IF( LCOEFFS ) ISEC4(6)  = 16
512C            `-----> wave fields have additional flags for
513C                    NC1 and NC2 coefficients (floats stored
514C                    in integer array isec4)
515      IRET = 1
516C
517C     If grid-point output, setup for 2nd order packing if requested.
518C
519      IF( (NOREPR.NE.JPSPHERE) .AND. (NOREPR.NE.JPSPHROT) ) THEN
520        HFUNC = 'C'
521        IF( NOHFUNC.EQ.'K' ) THEN
522          HFUNC = 'K'
523          ISEC4(4)  = 64
524          ISEC4(6)  = 16
525          ISEC4(9)  = 32
526          ISEC4(10) = 16
527          ISEC4(12) = 8
528          ISEC4(13) = 4
529          ISEC4(14) = 0
530          ISEC4(15) = -1
531        ELSE
532          ISEC4(4)  = 0
533        ENDIF
534      ELSE
535        HFUNC = 'C'
536        IF( NOHFUNC.EQ.'C' ) THEN
537          ISEC2(6) = 2
538          ISEC4(4) = 64
539        ELSE IF( NOHFUNC.EQ.'S' ) THEN
540          ISEC2(6) = 1
541          ISEC4(4) = 0
542        ENDIF
543
544      ENDIF
545
546C
547      CALL GRIBEX (ISEC0,ISEC1,ISEC2,ZSEC2,ISEC3,ZSEC3,ISEC4,
548     X             NEWAVE,ISEC4(1),OUTGRIB,OUTLEN,IWORD,HFUNC,IRET)
549C
550      IF ( IRET.NE.0 ) THEN
551        CALL INTLOG(JP_ERROR,'INTWAVE2: GRIBEX encoding failed.',IRET)
552        INTWAVE2 = JPROUTINE + 14
553        GOTO 900
554      ENDIF
555      OUTLEN = IWORD
556C
557C ------------------------------------------------------------------
558C*    Section 9.   Closedown.
559C ------------------------------------------------------------------
560C
561  900 CONTINUE
562C
563C     Clear change flags for next product processing
564C
565      LCHANGE = .FALSE.
566      LSMCHNG = .FALSE.
567C
568      RETURN
569      END
570