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 OCEANP( INOCEAN, INLEN, OTOCEAN, OUTLEN)
12C
13C---->
14C**** OCEANP
15C
16C     Purpose
17C     -------
18C
19C     Interpolate GRIB format input ocean field to GRIB format ocean
20C     field.
21C
22C
23C     Interface
24C     ---------
25C
26C     IRET = OCEANP( INOCEAN, INLEN, OTOCEAN, OUTLEN)
27C
28C     Input
29C     -----
30C
31C     INOCEAN - Input ocean field  (GRIB format).
32C     INLEN   - Input field length (words).
33C
34C
35C     Output
36C     ------
37C
38C     OTOCEAN - Output ocean field  (GRIB format).
39C     OUTLEN  - Output field length (words).
40C
41C
42C     Method
43C     ------
44C
45C     Interpolate ocean field.
46C
47C     If requested X grid increment does not match input GRIB
48C     increment, do a full interpolation in both x and y. Otherwise,
49C     interpolate only to make grid regular.
50C
51C     Externals
52C     ---------
53C
54C     JDEBUG  - Tests if debug output required.
55C     GRIBEX  - GRIB decoding/encoding.
56C     INTLOG  - Log error message.
57C     JMEMHAN - Allocate/deallocate scratch memory.
58C     INTOCN  - Carry out ocean field to field interpolation.
59C
60C
61C     Author
62C     ------
63C
64C     J.D.Chambers     ECMWF     May 1996
65C
66C     J.D.Chambers     ECMWF        Feb 1997
67C     Allow for 64-bit pointers
68C
69C     J.D.Chambers     ECMWF        Mar 1997
70C     Allow for full or partial interpolation
71C     (F or R option in call to INTOCN).
72C
73C----<
74C
75      IMPLICIT NONE
76C
77C     Function arguments
78C
79      INTEGER INOCEAN, INLEN, OTOCEAN, OUTLEN
80C
81#include "parim.h"
82#include "nifld.common"
83#include "nofld.common"
84#include "intf.h"
85C
86C     Parameters
87C
88      INTEGER JPROUTINE, JPALLOC, JPDEALL, JPSCR3, JPSCR4
89      INTEGER JP_LAT, JP_LONG, JP_GUESS
90      PARAMETER (JPROUTINE = 36900 )
91      PARAMETER (JPALLOC = 1)
92      PARAMETER (JPDEALL = 0)
93      PARAMETER (JPSCR3 = 3)
94      PARAMETER (JPSCR4 = 4)
95      PARAMETER (JP_LONG = 3)
96      PARAMETER (JP_LAT  = 4)
97      PARAMETER (JP_GUESS  = 18049691)
98cs    PARAMETER (JP_GUESS  = 1675245)
99cs    PARAMETER (JP_GUESS  = 1237104)
100cs    PARAMETER (JP_GUESS  = 1038240)
101C                          `---> allow for 0.25*0.25 interpolated grid
102C
103C     Local variables
104C
105      LOGICAL LHORIZN, LFLIP
106      CHARACTER*1 HFUNC
107      CHARACTER*1 HINTOPT
108      INTEGER IOLDN, IOLDW, IOLDS, IOLDE, IOLDEW, IOLDNS
109      INTEGER LOOPLAT, LOOPLON, I, J
110      REAL ZTEMP
111      INTEGER IERR, KPR, IWORD, IOLDSIZ, INEWSIZ, IINTPOL, ISTAGP
112      INTEGER NXP, NYP, INEWP, ITEMP
113      REAL RTOP, RBOTT, RRIGHT, RLEFT, RYRES, RXRES
114      REAL XARR, YARR
115      DIMENSION XARR(JPGRIB_ISEC1), YARR(JPGRIB_ISEC1)
116#ifndef _CRAYFTN
117#ifdef POINTER_64
118      INTEGER*8 INEW, IOLD
119#endif
120#endif
121      REAL NEW, OLD
122      POINTER ( INEW, NEW )
123      POINTER ( IOLD, OLD )
124      DIMENSION NEW( 1 ), OLD( 1 )
125C
126C     Externals
127C
128      INTEGER RESET_C, INTOCN, FIXAREA
129C
130C     -----------------------------------------------------------------|
131C*    Section 1.   Initialise
132C     -----------------------------------------------------------------|
133C
134  100 CONTINUE
135C
136C     Test if debug output required.
137C
138      CALL JDEBUG()
139C
140      CALL INTLOG(JP_DEBUG,
141     X  'OCEANP: Trying to interpolate an ocean field.',JPQUIET)
142C
143      IF( .NOT.LNOGRID ) THEN
144        CALL INTLOG(JP_FATAL,
145     X    'OCEANP: GRID must be specified to interpolate an ocean field'
146     X    , JPQUIET)
147        OCEANP = JPROUTINE + 1
148        GOTO 960
149      ENDIF
150C
151      OCEANP = 0
152      IERR   = 0
153      KPR    = 0
154      LFLIP = .FALSE.
155C
156C     Save the original area and grid definitions
157C
158      IOLDN  = NOAREA(1)
159      IOLDW  = NOAREA(2)
160      IOLDS  = NOAREA(3)
161      IOLDE  = NOAREA(4)
162      IOLDEW = NOGRID(1)
163      IOLDNS = NOGRID(2)
164C
165C     -----------------------------------------------------------------|
166C*    Section 2.   Unpack the input ocean field.
167C     -----------------------------------------------------------------|
168C
169  200 CONTINUE
170C
171C     Decode first four sections of GRIB code to find number of values
172C
173      IWORD   = INLEN
174      IERR    = 1
175      IOLDSIZ = OUTLEN
176      CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4,
177     X            OLD, IOLDSIZ, INOCEAN, INLEN, IWORD, 'I',IERR)
178C
179      IF( IERR.NE.0 ) THEN
180        CALL INTLOG(JP_ERROR,
181     X    'OCEANP: GRIBEX option J decoding failed.',IERR)
182        OCEANP = IERR
183        GOTO 950
184      ENDIF
185C
186C     Check it is an ocean field
187C
188      IF( (ISEC1(24).NE.1).OR.(ISEC1(37).NE.4) ) THEN
189        CALL INTLOG(JP_ERROR, 'OCEANP: Not an ocean field.',JPQUIET)
190        OCEANP = JPROUTINE + 3
191        GOTO 950
192      ENDIF
193C
194C     Horizontal (lat/long) field GRIB setup different from others
195C
196      LHORIZN = (ISEC1(60).EQ.JP_LONG).AND.(ISEC1(61).EQ.JP_LAT)
197      IF( LHORIZN ) CALL INTLOG(JP_DEBUG,
198     X  'OCEANP: Horizontal (lat/long) field interpolation.',JPQUIET)
199C
200C     Get scratch memory for unpacking the input ocean field.
201C
202      IOLDSIZ = ISEC2(2)*ISEC2(3)
203      CALL JMEMHAN( JPSCR4, IOLD, IOLDSIZ, JPALLOC, IERR)
204      IF( IERR.NE.0 ) THEN
205        CALL INTLOG(JP_ERROR,
206     X    'OCEANP: Scratch memory(4) allocation failed.',JPQUIET)
207        OCEANP = IERR
208        GOTO 950
209      ENDIF
210C
211C     Decode data from GRIB code: -999.9 is a missing data value
212C
213      IWORD = INLEN
214      IERR  = 1
215      ZSEC3(2) = -999.9
216      CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4,
217     X            OLD, IOLDSIZ, INOCEAN, INLEN, IWORD, 'D',IERR)
218C
219      IF( (IERR.NE.0).AND.(IERR.NE.-2).AND.(IERR.NE.-4) ) THEN
220        CALL INTLOG(JP_ERROR,'OCEANP: GRIBEX decoding failed.',IERR)
221        OCEANP = IERR
222        GOTO 940
223      ENDIF
224C
225C     -----------------------------------------------------------------|
226C*    Section 3.   Prepare the output ocean field.
227C     -----------------------------------------------------------------|
228C
229  300 CONTINUE
230C
231C     Setup interpolation options from input GRIB characteristics.
232C
233      IERR = RESET_C(ISEC1, ISEC2, ZSEC2, ISEC4)
234      IF( IERR.NE.0 ) THEN
235        CALL INTLOG(JP_ERROR,
236     X    'OCEANP: Setup interp. options from GRIB failed.',JPQUIET)
237        OCEANP = IERR
238        GOTO 940
239      ENDIF
240C
241C     For horizontal field, match the desired area to the grid
242C
243      IERR = FIXAREA()
244      IF( IERR.NE.0 ) THEN
245        CALL INTLOG(JP_ERROR, 'OCEANP: FIXAREA failed.',JPQUIET)
246        OCEANP = IERR
247        GOTO 940
248      ENDIF
249C
250C     Get scratch memory for output ocean field.
251C     Have to guess how big the interpolated grid will be.
252C
253      INEWSIZ = JP_GUESS
254      CALL JMEMHAN(JPSCR3, INEW, INEWSIZ*2, JPALLOC, IERR)
255      IF( IERR.NE.0 ) THEN
256        CALL INTLOG(JP_ERROR,
257     X    'OCEANP: Scratch memory(3) allocation failed.',JPQUIET)
258        OCEANP = IERR
259        GOTO 940
260      ENDIF
261      INEWP = 1
262      ITEMP = 1 + INEWSIZ
263C
264C     -----------------------------------------------------------------|
265C*    Section 4.   Interpolate ocean field.
266C     -----------------------------------------------------------------|
267C
268  400 CONTINUE
269C
270      RTOP   = FLOAT( NOAREA(1) ) / PPMULT
271      RBOTT  = FLOAT( NOAREA(3) ) / PPMULT
272      RRIGHT = FLOAT( NOAREA(4) ) / PPMULT
273      RLEFT  = FLOAT( NOAREA(2) ) / PPMULT
274C
275      RYRES = FLOAT( NOGRID(2) ) / PPMULT
276      RXRES = FLOAT( NOGRID(1) ) / PPMULT
277C
278C     Select the interpolation option
279C     (default = 'R' = interpolate only to make grid regular)
280C     If requested X grid increment does not match input GRIB
281C     increment, interpolate in x and y ('F' = full)
282C
283      HINTOPT = 'R'
284      IF( NOGRID(1).NE.ISEC2(9) ) HINTOPT = 'F'
285C
286      IERR = INTOCN(HINTOPT,'D',
287     X              RTOP, RBOTT, RRIGHT, RLEFT, RYRES, RXRES,
288     X              ISEC1, ISEC2, OLD, IOLDSIZ,
289     X              INEWSIZ, NEW(ITEMP), XARR, YARR, NXP, NYP,
290     X              NEW(INEWP),
291     X              IINTPOL, ISTAGP)
292      IF( IERR.NE.0 ) THEN
293        CALL INTLOG(JP_ERROR, 'OCEANP: Interpolation failed.',JPQUIET)
294        OCEANP = IERR
295        GOTO 900
296      ENDIF
297C
298      NOAREA(1) = NINT( (YARR(1) + YARR(2)*(NYP-1)) * PPMULT )
299      NOAREA(2) = NINT( XARR(1) * PPMULT )
300      NOAREA(3) = NINT( YARR(1) * PPMULT )
301      NOAREA(4) = NINT( (XARR(1) + XARR(2)*(NXP-1)) * PPMULT )
302C
303      NOGRID(1) = NINT( ABS(XARR(2)) * PPMULT )
304      NOGRID(2) = NINT( ABS(YARR(2)) * PPMULT )
305      NONS = NYP
306      NOWE = NXP
307C
308C     -----------------------------------------------------------------|
309C*    Section 5.   Code new field into GRIB.
310C     -----------------------------------------------------------------|
311C
312  500 CONTINUE
313C
314C     Re-order the rows from south to north to be north to south
315C     (ECMWF convention).
316C
317      IF( LHORIZN.OR.(NOAREA(1).LT.NOAREA(3)) ) THEN
318        LFLIP = .TRUE.
319        DO LOOPLAT = 1, NYP/2
320          I = (LOOPLAT - 1)*NXP
321          J = (NYP - LOOPLAT)*NXP
322          DO LOOPLON = 1, NXP
323            I = I + 1
324            J = J + 1
325            ZTEMP  = NEW(I)
326            NEW(I) = NEW(J)
327            NEW(J) = ZTEMP
328          ENDDO
329        ENDDO
330C
331        IF( .NOT.LHORIZN ) THEN
332          ITEMP     = NOAREA(1)
333          NOAREA(1) = NOAREA(3)
334          NOAREA(3) = ITEMP
335C
336          ITEMP    = ISEC2(4)
337          ISEC2(4) = ISEC2(7)
338          ISEC2(7) = ITEMP
339        ENDIF
340      ENDIF
341C
342C     Setup GRIB sections for the output product:
343C
344C     Use local definition 4 in GRIB section 1 ...
345C
346      ISEC1(73) = 0
347      ITEMP = 75 + ISEC1(71) + ISEC1(72) + ISEC1(73) + ISEC1(74)
348      ISEC1(ITEMP) = 0
349
350C
351      IF( LHORIZN ) THEN
352        ISEC1(62) = NOAREA(1)*10
353        ISEC1(63) = NOAREA(2)*10
354        ISEC1(64) = NOAREA(3)*10
355        ISEC1(65) = NOAREA(4)*10
356        ISEC1(66) = NOGRID(1)*10
357        ISEC1(67) = NOGRID(2)*10
358        IF( LFLIP ) ISEC1(67) = - ISEC1(67)
359        ISEC1(68) = 0
360        ISEC1(69) = 0
361C
362C       Standard section 2 format for horizontal sections
363C
364        ISEC2(1) = 0
365        ISEC2(2) = NXP
366        ISEC2(3) = NYP
367        ISEC2(4) = NOAREA(1)/100
368        ISEC2(5) = NOAREA(2)/100
369        ISEC2(6) = 128
370        ISEC2(7) = NOAREA(3)/100
371        ISEC2(8) = NOAREA(4)/100
372        ISEC2(9)  = NOGRID(1)/100
373        ISEC2(10) = NOGRID(2)/100
374Cjdc    ISEC2(17) = 0
375        DO LOOPLAT = 11, 22
376          ISEC2(LOOPLAT) = 0
377        ENDDO
378      ELSE
379C Tim-Correct headers for vertical sections
380        IF( ISEC1(60).EQ.3 ) THEN
381          ISEC1(63) = NOAREA(2)*10
382        ELSEIF( ISEC1(60).EQ.4 ) THEN
383          ISEC1(63) = NOAREA(2)*10
384        ENDIF
385        IF( ISEC1(61).EQ.3 ) THEN
386          ISEC1(62) = NOAREA(1)*10
387        ELSEIF( ISEC1(61).EQ.4 ) THEN
388          ISEC1(62) = NOAREA(1)*10
389        ENDIF
390C
391        IF( ISEC1(61).EQ.1 ) THEN
392          ISEC1(64) = ISEC1(62) + NINT((NYP-1)*RYRES*1.0)
393          ISEC1(65) = ISEC1(63) + NINT((NXP-1)*RXRES*1000000.0)
394          ISEC1(66) = NOGRID(1)*NINT(1000000.0/PPMULT)
395          ISEC1(67) = NOGRID(2)/(PPMULT/1.0)
396        ELSE
397          ISEC1(64) = ISEC1(62) + NINT((NYP-1)*RYRES*1000.0)
398          ISEC1(65) = ISEC1(63) + NINT((NXP-1)*RXRES*1000000.0)
399          ISEC1(66) = NOGRID(1)*NINT(1000000.0/PPMULT)
400          ISEC1(67) = NOGRID(2)/(PPMULT/1000.0)
401        ENDIF
402C
403        IF( LFLIP ) ISEC1(67) = - ISEC1(67)
404        ISEC1(68) = 0
405C
406C       Special ECMWF ocean short form of GRIB section 2
407C
408        ISEC2(1) = 192
409        ISEC2(2) = NXP
410        ISEC2(3) = NYP
411Cjdc    ISEC2(4) = NOAREA(1)/100000
412Cjdc    ISEC2(5) = NOAREA(2)/100
413Cjdc    ISEC2(6) = 0
414Cjdc    ISEC2(7) = NOAREA(3)/100000
415Cjdc    ISEC2(8) = NOAREA(4)/100
416Cjdc    DO LOOPLAT = 9, 22
417Cjdc      ISEC2(LOOPLAT) = 0
418Cjdc    ENDDO
419      ENDIF
420C
421C     Ensure there is a section 3 in case there is missing data in the
422C     output area
423C
424      ISEC1(5) = 192
425      ISEC3(1) = 0
426      ZSEC3(2) = -999.9
427C
428      ISEC4(1) = NXP *  NYP
429C
430C     If grid-point output, setup for 2nd order packing if requested.
431C
432      IF( (NOREPR.NE.JPSPHERE).AND.(NOREPR.NE.JPSPHROT) ) THEN
433        HFUNC = 'C'
434        IF( NOHFUNC.EQ.'K' ) THEN
435          HFUNC = 'K'
436          ISEC4(4)  = 64
437          ISEC4(6)  = 16
438          ISEC4(9)  = 32
439          ISEC4(10) = 16
440          ISEC4(12) = 8
441          ISEC4(13) = 4
442          ISEC4(14) = 0
443          ISEC4(15) = -1
444        ELSE
445          ISEC4(4)  = 0
446          ISEC4(6)  = 0
447        ENDIF
448      ELSE
449        HFUNC = 'C'
450        IF( NOHFUNC.EQ.'C' ) THEN
451          ISEC2(6) = 2
452          ISEC4(4) = 64
453        ELSE IF( NOHFUNC.EQ.'S' ) THEN
454          ISEC2(6) = 1
455          ISEC4(4) = 0
456        ENDIF
457      ENDIF
458C
459C     Switch off GRIBEX checking (sigh)
460C
461      CALL INTLOG(JP_DEBUG,
462     X  'OCEANP: switching OFF GRIBEX checking',JPQUIET)
463      CALL GRSVCK(0)
464C
465      IERR = 1
466      CALL GRIBEX( ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4,
467     X             NEW,INEWSIZ,OTOCEAN,OUTLEN,IWORD,HFUNC,IERR)
468C
469      IF( IERR.NE.0 ) THEN
470        CALL INTLOG(JP_ERROR,'OCEANP: GRIBEX encoding failed.',IERR)
471        OCEANP = JPROUTINE + 5
472        GOTO 900
473      ENDIF
474      OUTLEN = IWORD
475C
476C     Switch on GRIBEX checking (sigh again)
477C
478      CALL INTLOG(JP_DEBUG,
479     X  'OCEANP: switching ON GRIBEX checking',JPQUIET)
480      CALL GRSVCK(1)
481C
482C     -----------------------------------------------------------------|
483C*    Section 9.   Closedown.
484C     -----------------------------------------------------------------|
485C
486  900 CONTINUE
487C
488C     Clear change flags for next product processing
489C
490      LCHANGE = .FALSE.
491      LSMCHNG = .FALSE.
492C
493C     Return the scratch memory.
494C
495      CALL JMEMHAN( JPSCR3, INEW, INEWSIZ, JPDEALL, IERR)
496      IF ( IERR .NE. 0 ) THEN
497        CALL INTLOG(JP_ERROR,
498     X    'OCEANP: Scratch memory(3) reallocation failed.',JPQUIET)
499        OCEANP = IERR
500      ENDIF
501C
502  940 CONTINUE
503C
504      CALL JMEMHAN( JPSCR4, IOLD, IOLDSIZ, JPDEALL, IERR)
505      IF ( IERR .NE. 0 ) THEN
506        CALL INTLOG(JP_ERROR,
507     X    'OCEANP: Scratch memory(4) reallocation failed.',JPQUIET)
508        OCEANP = IERR
509      ENDIF
510C
511  950 CONTINUE
512C
513C     Restore the original area and grid definitions
514C
515      NOAREA(1) = IOLDN
516      NOAREA(2) = IOLDW
517      NOAREA(3) = IOLDS
518      NOAREA(4) = IOLDE
519      NOGRID(1) = IOLDEW
520      NOGRID(2) = IOLDNS
521C
522  960 CONTINUE
523C
524      CALL INTLOG(JP_DEBUG,
525     X  'OCEANP: Returning from interpolating an ocean field.',JPQUIET)
526C
527      RETURN
528      END
529