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 HNTFAPH(INGRIB,INLEN)
12C
13C---->
14C**** HNTFAPH
15C
16C     Purpose
17C     -------
18C
19C     Prepare to interpolate to grid point field.
20C
21C
22C     Interface
23C     ---------
24C
25C     IRET = HNTFAPH(INGRIB,INLEN)
26C
27C     Input
28C     -----
29C
30C     INGRIB - Input field (packed).
31C     INLEN  - Input field length (words).
32C
33C
34C     Output
35C     ------
36C
37C     Field unpacked values are in ZNFELDI, rotated if necessary.
38C
39C     Returns: 0, if OK. Otherwise, an error occured in interpolation.
40C
41C
42C     Method
43C     ------
44C
45C     Unpack field.
46C
47C     If the input is a spectral field and the output is a rotated
48C     grid-point field, create a global reduced gaussian field and
49C     then create the rotated grid-point field from it.
50C
51C
52C     Externals
53C     ---------
54C
55C     GRIBEX  - Decode/encode GRIB product.
56C     GRSVCK  - Turn off GRIB checking
57C     INTLOG  - Log error message.
58C     INTLOGR - Log error message.
59C     JDEBUG  - Checks environment to switch on/off debug
60C     FIXAREA - Fixup input/output field area definitions.
61C     HSH2GG  - Find suitable gaussian grid/spectral truncation
62C     HIRLAM  - Creates rotated lat/long field from reduced gaussian
63C     HIRLSM  - Creates rotated lat/long field from reduced gaussian
64C               using land-sea mask
65C     HRG2GG  - Creates rotated gaussian field from reduced gaussian
66C     HLL2LL  - Creates rotated lat/long field from lat/long field
67C     LSMFLD  - Determines whether a field is to be interpolated using
68C               a land-sea mask
69C
70C
71C     Author
72C     ------
73C
74C     J.D.Chambers     ECMWF     January 31, 2001
75C
76C----<
77C
78C     -----------------------------------------------------------------|
79C*    Section 0.   Variables
80C     -----------------------------------------------------------------|
81C
82      IMPLICIT NONE
83C
84C     Function arguments
85C
86      INTEGER INGRIB(*),INLEN
87C
88#include "parim.h"
89#include "nifld.common"
90#include "nofld.common"
91#include "intf.h"
92#include "current.h"
93#include "grfixed.h"
94C
95C     Parameters
96C
97      INTEGER JPROUTINE
98      PARAMETER (JPROUTINE = 40110 )
99C
100C     Local variables
101C
102      CHARACTER*1 HOLDTYP
103      CHARACTER*1 HTYPE
104      INTEGER     IERR
105      INTEGER     IRET
106      INTEGER     ISHIZE
107      INTEGER     ISIZE
108      INTEGER     ITEMP,K,I,J,IM,JM
109      INTEGER     IWORD
110      INTEGER     KPTS(JPGTRUNC*2)
111      INTEGER     LOOP
112      INTEGER     NCOUNT
113      INTEGER     NGAUSS
114      INTEGER     NLAT
115      INTEGER     NLON
116      INTEGER     NTRUNC
117      INTEGER     NUMPTS
118      INTEGER     NUVFLAG
119      LOGICAL     LLATOUT
120      LOGICAL     LSP2RGG
121      LOGICAL     LUSELSM
122      REAL        AREA(4)
123      REAL        EAST
124      REAL        EW
125      REAL        GRID(2)
126      REAL        NS
127      REAL        OLDGRID(2)
128      REAL        PLATS(JPGTRUNC*2)
129      REAL        POLE(2)
130      REAL        RGGRID(1)
131      REAL        SWORK(1)
132      REAL        TEMP(1440,1440)
133      REAL        WEST
134      REAL        ZNFLDO(1)
135#ifndef _CRAYFTN
136#ifdef POINTER_64
137      INTEGER*8 IRGGRID
138      INTEGER*8 ISWORK
139      INTEGER*8 IZNFLDO
140#endif
141#endif
142      POINTER (IRGGRID, RGGRID)
143      POINTER (ISWORK,  SWORK)
144      POINTER (IZNFLDO, ZNFLDO)
145C
146      SAVE IRGGRID
147      SAVE ISWORK
148      SAVE IZNFLDO
149C
150C     Externals
151C
152      CHARACTER*1 GGHTYPE
153      INTEGER HSH2GG, HIRLAM, HIRLSM, HRG2GG, FIXAREA, PDDEFS
154      INTEGER HLL2LL
155      LOGICAL LSMFLD
156      EXTERNAL GGHTYPE
157      EXTERNAL HSH2GG, HIRLAM, HIRLSM, HRG2GG, FIXAREA, PDDEFS
158      EXTERNAL HLL2LL
159      EXTERNAL LSMFLD
160C
161C     -----------------------------------------------------------------|
162C*    Section 1.   Initialise
163C     -----------------------------------------------------------------|
164C
165  100 CONTINUE
166C
167      HNTFAPH = 0
168      LSP2RGG = .FALSE.
169C
170C     -----------------------------------------------------------------|
171C*    Section 2.   Decode data from the GRIB code
172C     -----------------------------------------------------------------|
173C
174  200 CONTINUE
175C
176C     Decode data from GRIB code (no checking)
177C
178      IERR  =  0
179      CALL GRSVCK(0)
180C
181      IERR = 1
182      IWORD = INLEN
183      ISEC3(2) = NINT(RMISSGV)
184      ZSEC3(2) = RMISSGV
185
186      CALL GRIBEX(ISEC0, ISEC1, ISEC2, ZSEC2, ISEC3, ZSEC3, ISEC4,
187     X            ZNFELDI, JPEXPAND, INGRIB, INLEN, IWORD, 'D',IERR)
188
189      IF (ISEC2(11).EQ.64) THEN
190        CALL INTLOG(JP_DEBUG,
191     X    'HNTFAPH: Scanning flag west-east/south-north',ISEC2(11))
192        ITEMP = NIAREA(1)
193        NIAREA(1) = NIAREA(3)
194        NIAREA(3) = ITEMP
195
196        IM = ISEC2(2)
197        JM = ISEC2(3)
198        K=0
199        DO J=JM,1,-1
200          DO I=1,IM
201            K=K+1
202            TEMP(I,J) = ZNFELDI(K)
203          END DO
204        END DO
205        K=0
206        DO J=1,JM
207          DO I=1,IM
208            K=K+1
209            ZNFELDI(K) = TEMP(I,J)
210          END DO
211        END DO
212
213      ENDIF
214
215C
216      IF( IERR.LT.0) THEN
217        IF( (IERR.EQ.-2).OR.(IERR.EQ.-4) ) THEN
218          CALL INTLOG(JP_DEBUG,'HNTFAPH: Use missing value',JPQUIET)
219          LIMISSV = .TRUE.
220        ELSE
221          CALL INTLOG(JP_ERROR,'HNTFAPH: GRIBEX decoding fail.',IERR)
222          HNTFAPH = JPROUTINE + 2
223          GOTO 900
224        ENDIF
225      ELSE IF( IERR.GT.0 ) THEN
226        CALL INTLOG(JP_ERROR,'HNTFAPH: GRIBEX decoding failed.',IERR)
227        HNTFAPH = JPROUTINE + 2
228        GOTO 900
229      ENDIF
230C
231      NCOUNT = ISEC4(1)
232C
233      LLATOUT = (NOREPR.EQ.JPREGROT).OR.(NOREPR.EQ.JPREGULAR)
234
235      IF( .NOT.LNOROTA ) GOTO 900
236C
237C     -----------------------------------------------------------------|
238C*    Section 3.   Handle rotation, if necessary.
239C     -----------------------------------------------------------------|
240C
241  300 CONTINUE
242C
243      CALL INTLOG(JP_DEBUG,'HNTFAPH: Rotate field.',JPQUIET)
244      CALL INTLOG(JP_DEBUG,'HNTFAPH: South pole lat  ',NOROTA(1))
245      CALL INTLOG(JP_DEBUG,'HNTFAPH: South pole long ',NOROTA(2))
246C
247C     Fill area limits (handles case when default 0/0/0/0 given)
248C
249      IRET = FIXAREA()
250      IF( IRET.NE.0 ) THEN
251        CALL INTLOG(JP_ERROR,'HNTFAPH: area fixup failed',JPQUIET)
252        HNTFAPH = JPROUTINE + 3
253        GOTO 900
254      ENDIF
255C
256      AREA(1) = REAL(NOAREA(1))/PPMULT
257      AREA(2) = REAL(NOAREA(2))/PPMULT
258      AREA(3) = REAL(NOAREA(3))/PPMULT
259      AREA(4) = REAL(NOAREA(4))/PPMULT
260C
261      GRID(1) = REAL(NOGRID(1))/PPMULT
262      GRID(2) = REAL(NOGRID(2))/PPMULT
263C
264      POLE(1) = REAL(NOROTA(1))/PPMULT
265      POLE(2) = REAL(NOROTA(2))/PPMULT
266C
267C     -----------------------------------------------------------------|
268C*    Section 4.   Spectral to rotated grid-point
269C     -----------------------------------------------------------------|
270C
271  400 CONTINUE
272C
273      IF( (NIREPR.EQ.JPSPHERE).OR.(NIREPR.EQ.JPSPHROT) ) THEN
274C
275C       Convert spectral to suitable global reduced gaussian grid
276C
277        CALL INTLOG(JP_DEBUG,
278     X    'HNTFAPH: Spectral to suitable reduced gaussian',JPQUIET)
279C
280        NIRESO = ISEC2(2)
281        NTRUNC = ISEC2(2)
282        IF( LNORESO ) THEN
283          NTRUNC = NORESO
284          NGAUSS = 0
285          HTYPE  = ''
286          NS = 0.
287          EW = 0.
288          IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE)
289          IF( IRET.NE.0 ) THEN
290            CALL INTLOG(JP_ERROR,
291     X        'HNTFAPH: problem getting data for reduced grid',NTRUNC)
292            HNTFAPH = JPROUTINE + 4
293            GOTO 900
294          ENDIF
295          GOTO 401
296        ENDIF
297
298        IF( LARESOL ) THEN
299          NS = 0.
300          EW = 0.
301          IF( LLATOUT ) THEN
302            NS = GRID(1)
303            EW = GRID(2)
304            NTRUNC = 0
305            NGAUSS = 0
306            HTYPE = ''
307          ELSE
308            HTYPE = 'R'
309          ENDIF
310          IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE)
311          IF( IRET.NE.0 ) THEN
312            CALL INTLOG(JP_ERROR,
313     X        'HNTFAPH: problem getting data for reduced grid',NTRUNC)
314            HNTFAPH = JPROUTINE + 4
315            GOTO 900
316          ENDIF
317        ENDIF
318
319C       Truncate if a smaller resolution has been requested
320C
321  401 CONTINUE
322      IF( NTRUNC.LT.NIRESO ) THEN
323       CALL INTLOG(JP_DEBUG,'HNTFAPH: Truncation changed from:',NIRESO)
324       CALL INTLOG(JP_DEBUG,'HNTFAPH: to: ',NTRUNC)
325       CALL INTLOG(JP_DEBUG,'HNTFAPH: Gaussian number is:',NGAUSS)
326C
327        ISHIZE =  (NTRUNC+1)*(NTRUNC+4)
328        CALL JMEMHAN( 3, IZNFLDO, ISHIZE, 1, IERR)
329        IF( IERR.NE.0 ) THEN
330          CALL INTLOG(JP_FATAL,
331     X      'HNTFAPH: Get scratch space failed',JPQUIET)
332          HNTFAPH = JPROUTINE + 4
333          GOTO 900
334        ENDIF
335
336C       Generate spherical harmonics with output truncation
337        CALL SH2SH( ZNFELDI, NIRESO, ZNFLDO, NTRUNC )
338
339C       Move new spherical harmonics to 'input' array
340        ZNFELDI(1:ISHIZE) = ZNFLDO(1:ISHIZE)
341
342      ELSE
343          CALL INTLOG(JP_DEBUG,
344     X      'HNTFAPH: Spectral to suitable reduced gaussian',JPQUIET)
345C
346          NTRUNC = ISEC2(2)
347          NGAUSS = 0
348          HTYPE  = ''
349          NS = 0.
350          EW = 0.
351          IRET = HSH2GG(NS,EW,NTRUNC,NGAUSS,HTYPE,KPTS,PLATS,ISIZE)
352          IF( IRET.NE.0 ) THEN
353            CALL INTLOG(JP_ERROR,
354     X        'HNTFAPH: problem getting data for reduced grid',NTRUNC)
355            HNTFAPH = JPROUTINE + 4
356            GOTO 900
357          ENDIF
358
359        ENDIF
360
361C
362C       Dynamically allocate memory for global reduced gaussian grid
363C
364        CALL JMEMHAN( 18, IRGGRID, ISIZE, 1, IRET)
365        IF( IRET.NE.0 ) THEN
366          CALL INTLOG(JP_ERROR,
367     X      'HNTFAPH: memory alloc for reduced grid fail',JPQUIET)
368          HNTFAPH = JPROUTINE + 4
369          GOTO 900
370        ENDIF
371
372C       Set flag to show field is not a wind component
373        NUVFLAG = 0
374
375C       Create the reduced gaussian grid
376        HOLDTYP = HOGAUST
377        WEST = 0.0
378        EAST = 360.0 - (360.0/(NGAUSS*4))
379c       EMOS-199: adjusted for reduced_gg/octahedral
380        IF (HTYPE.EQ.'O')  EAST = 360.0 - (360.0/FLOAT(KPTS(NGAUSS)))
381        IF (HTYPE.NE.'R' .AND. HTYPE.NE.'O' .AND.
382     X      HTYPE.NE.'F' .AND. HTYPE.NE.'U') THEN
383          HTYPE = 'R'
384        ENDIF
385        CALL JAGGGP(ZNFELDI,NTRUNC,PLATS(1),PLATS(NGAUSS*2),WEST,
386     X              EAST,NGAUSS,HTYPE,KPTS,RGGRID,NUVFLAG,IRET)
387        IF( IRET.NE.0 ) THEN
388          CALL INTLOG(JP_ERROR,
389     X      'HNTFAPH: spectral to reduced gaussian failed',JPQUIET)
390          HNTFAPH = JPROUTINE + 4
391          GOTO 900
392        ENDIF
393        HOGAUST = HOLDTYP
394
395        NCOUNT = 0
396        DO LOOP = 1, (NGAUSS*2)
397          NCOUNT= NCOUNT + KPTS(LOOP)
398        ENDDO
399
400        LSP2RGG = .TRUE.
401
402        IF( NOREPR.EQ.JPFGGROT ) THEN
403          CALL INTLOG(JP_DEBUG,
404     X      'HNTFAPH: Convert gaussian to rotated gaussian',JPQUIET)
405          GOTO 600
406        ENDIF
407
408      ENDIF
409C
410C     -----------------------------------------------------------------|
411C*    Section 5.   Gaussian to rotated lat/long
412C     -----------------------------------------------------------------|
413C
414  500 CONTINUE
415C
416      IF( (LLATOUT) .AND. (
417     X        (NIREPR.EQ.JPQUASI) .OR.
418     X        (NIREPR.EQ.JPGAUSSIAN) .OR.
419     X        (LSP2RGG) )) THEN
420C
421        CALL INTLOG(JP_DEBUG,'HNTFAPH: Gauss to lat/lon',JPQUIET)
422C
423C       Dynamically allocate memory for rotated lat/long grid
424C
425        NLON = 1 + NINT((AREA(JPEAST)  - AREA(JPWEST)) /
426     X         GRID(JPWESTEP)) !SC
427        NLAT = 1 + NINT((AREA(JPNORTH) - AREA(JPSOUTH)) /
428     X         GRID(JPNSSTEP)) !SC
429C
430        NOWE = NLON
431        NONS = NLAT
432C
433        NUMPTS = NLON * NLAT
434        ISIZE  = NUMPTS
435        CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET)
436        IF( IRET.NE.0 ) THEN
437          CALL INTLOG(JP_ERROR,
438     X      'HNTFAPH: memory alloc for lat/long grid fail',JPQUIET)
439          HNTFAPH = JPROUTINE + 5
440          GOTO 900
441        ENDIF
442C
443        LUSELSM = LSMFLD()
444C
445C       If original field was spectral, ...
446C
447        IF( LSP2RGG ) THEN
448          IF( LUSELSM ) THEN
449            IRET = HIRLSM(LO12PT,RGGRID,NCOUNT,NGAUSS,HTYPE,AREA,POLE,
450     X                    GRID,SWORK,ISIZE,NLON,NLAT)
451          ELSE
452            IRET = HIRLAM(LO12PT,RGGRID,NCOUNT,NGAUSS,HTYPE,AREA,POLE,
453     X                    GRID,SWORK,ISIZE,NLON,NLAT)
454          ENDIF
455C
456        ELSE
457C
458C       If original field was gaussian, ...
459C
460          IRET = PDDEFS()
461          NGAUSS = ISEC2(10)
462          HTYPE  = GGHTYPE(NIREPR,NGAUSS,MILLEN)
463cs        IF( LUSELSM ) THEN
464          IF( LSM ) THEN
465            IRET = HIRLSM(LO12PT,ZNFELDI,NCOUNT,NGAUSS,HTYPE,AREA,POLE,
466     X                    GRID,SWORK,ISIZE,NLON,NLAT)
467          ELSE
468            IRET = HIRLAM(LO12PT,ZNFELDI,NCOUNT,NGAUSS,HTYPE,AREA,POLE,
469     X                    GRID,SWORK,ISIZE,NLON,NLAT)
470          ENDIF
471C
472        ENDIF
473C
474        IF( IRET.NE.0 ) THEN
475          CALL INTLOG(JP_ERROR,
476     X      'HNTFAPH: HIRLAM rotation failed',JPQUIET)
477          HNTFAPH = JPROUTINE + 7
478          GOTO 900
479        ENDIF
480C
481        ISEC2(1) = JPREGROT
482        ISEC4(1) = NOWE * NONS
483C
484      ENDIF
485C
486C     -----------------------------------------------------------------|
487C*    Section 6.   Gaussian to rotated gaussian
488C     -----------------------------------------------------------------|
489C
490  600 CONTINUE
491C
492      IF( (LSP2RGG.AND.(NOREPR.EQ.JPFGGROT)) .OR.
493     X    (LSP2RGG.AND.(NOREPR.EQ.JPQGGROT)) .OR.
494     X    (((NIREPR.EQ.JPQUASI).OR.(NIREPR.EQ.JPGAUSSIAN)) .AND.
495     X     (NOREPR.EQ.JPFGGROT).OR.(NOREPR.EQ.JPQGGROT)) ) THEN
496        CALL INTLOG(JP_DEBUG,
497     X    'HNTFAPH: Gaussian to reduced gaussian',JPQUIET)
498C
499C       Dynamically allocate memory for rotated lat/long grid
500C
501C       ISIZE = NOGAUSS * NOGAUSS * 8
502        ISIZE = (2*NOGAUSS) * (4*NOGAUSS + 20)  ! account for RGG/octahedral
503        CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET)
504        IF( IRET.NE.0 ) THEN
505          CALL INTLOG(JP_ERROR,
506     X      'HNTFAPH: memory alloc for gaussian grid fail',JPQUIET)
507          HNTFAPH = JPROUTINE + 6
508          GOTO 900
509        ENDIF
510C
511C       If original field was spectral, ...
512C
513        IF( LSP2RGG ) THEN
514          IRET = HRG2GG(LO12PT,RGGRID,NGAUSS,AREA,POLE,
515     X                  NOGAUSS,HOGAUST,SWORK,ISIZE,NUMPTS)
516C
517        ELSE
518C
519C       If original field was gaussian, ...
520C
521          NGAUSS = ISEC2(10)
522          IRET = HRG2GG(LO12PT,ZNFELDI,NGAUSS,AREA,POLE,
523     X                  NOGAUSS,HOGAUST,SWORK,ISIZE,NUMPTS)
524        ENDIF
525        IF( IRET.NE.0 ) THEN
526          CALL INTLOG(JP_ERROR,
527     X      'HNTFAPH: HRG2GG rotation failed',JPQUIET)
528          HNTFAPH = JPROUTINE + 6
529          GOTO 900
530        ENDIF
531C
532        IF( (NOREPR.EQ.JPQUASI).OR.
533     X      (NOREPR.EQ.JPQGGROT).OR.
534     X      (NOREPR.EQ.JPFGGROT).OR.
535     X      (NOREPR.EQ.JPGAUSSIAN) ) THEN
536          ISEC2(1) = JPFGGROT
537        ELSE
538          ISEC2(1) = NOREPR
539        ENDIF
540C
541        ISEC4(1) = NUMPTS
542C
543      ENDIF
544C
545C     -----------------------------------------------------------------|
546C*    Section 7.   Lat/long to rotated lat/long
547C     -----------------------------------------------------------------|
548C
549  700 CONTINUE
550C
551      IF( (NIREPR.EQ.JPREGULAR) ) THEN
552C
553C       Dynamically allocate memory for rotated lat/long grid
554C
555        NLON = 1 + NINT((AREA(JPEAST)  - AREA(JPWEST)) /
556     X         GRID(JPWESTEP))
557        NLAT = 1 + NINT((AREA(JPNORTH) - AREA(JPSOUTH)) /
558     X         GRID(JPNSSTEP))
559C
560        NOWE = NLON
561        NONS = NLAT
562C
563        NUMPTS = NLON * NLAT
564        ISIZE  = NUMPTS
565        CALL JMEMHAN( 11, ISWORK, ISIZE, 1, IRET)
566        IF( IRET.NE.0 ) THEN
567          CALL INTLOG(JP_ERROR,
568     X      'HNTFAPH: memory alloc for lat/long grid fail',JPQUIET)
569          HNTFAPH = JPROUTINE + 7
570          GOTO 900
571        ENDIF
572C
573        OLDGRID(1) = REAL(NIGRID(1))/100000.0
574        OLDGRID(2) = REAL(NIGRID(2))/100000.0
575        IRET = HLL2LL(LO12PT,ZNFELDI,OLDGRID,AREA,POLE,GRID,SWORK,ISIZE,
576     X                NLON,NLAT)
577C
578        IF( IRET.NE.0 ) THEN
579          CALL INTLOG(JP_ERROR,
580     X      'HNTFAPH: HLL2LL rotation failed',JPQUIET)
581          HNTFAPH = JPROUTINE + 7
582          GOTO 900
583        ENDIF
584C
585        ISEC2(1) = JPREGROT
586        ISEC4(1) = NOWE * NONS
587C
588      ENDIF
589
590C     move rotated field back into field original array
591      ZNFELDI(1:NUMPTS) = SWORK(1:NUMPTS)
592
593C     -----------------------------------------------------------------|
594C*    Section 9.   Closedown.
595C     -----------------------------------------------------------------|
596C
597  900 CONTINUE
598C
599      RETURN
600      END
601