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 HRG2GG(L12PNT,OLDFLD,KGNOLD,AREA,POLE,
12     X                        KGNNEW,HTYPE,NEWFLD,KSIZE,NUMPTS)
13C
14C---->
15C**** HRG2GG
16C
17C     Purpose
18C     -------
19C
20C     This routine creates a rotated gaussian field from a reduced
21C     gaussian field using 12-point horizontal interpolation.
22C
23C
24C     Interface
25C     ---------
26C
27C     IRET = HRG2GG(L12PNT,OLDFLD,KGNOLD,AREA,POLE,
28C    X              KGNNEW,HTYPE,NEWFLD,KSIZE,NUMPTS)
29C
30C
31C     Input parameters
32C     ----------------
33C
34C     L12PNT  - Chooses between 12-point and 4-point interpolation
35C     OLDFLD  - Array of values from the old (reduced) gaussian field
36C     KGNOLD  - Gaussian number for the old gaussian field
37C     AREA    - Limits of area (N/W/S/E)
38C     POLE    - Pole of rotation (lat/long)
39C     KGNNEW  - Gaussian grid number (for the new field)
40C     HTYPE   - Gaussian grid type (for the new field)
41C               = 'R' for reduced ("quasi-regular"),
42C               = 'O' for reduced/octahedral,
43C               = 'F' for full
44C     KSIZE   - The size of the array to fill with the regular
45C               lat/long field
46C
47C
48C     Output parameters
49C     -----------------
50C
51C     NEWFLD  - The array of values for the gaussian field
52C     NUMPTS  - Number of points in the new gaussian field.
53C
54C     Returns 0 if function successful, non-zero otherwise.
55C
56C
57C     Common block usage
58C     ------------------
59C
60C     None
61C
62C
63C     Method
64C     ------
65C
66C     Numbering of the points (I is the interpolation point):
67C
68C                   13       5       6      14
69C
70C                    7       1       2       8
71C                               (I)
72C                    9       3       4      10
73C
74C                   15      11      12      16
75C
76C     The 12-point interpolation is not possible if either of the top
77C     two rows is above the original field northern latitude. The
78C     nearest neighbour is used if both rows are above, and a 4-pt
79C     bilinear interpolation is used if the top row is above.
80C     Similarily, if either of the bottom two rows is below the original
81C     field southern latitude.
82C
83C
84C     Externals
85C     ---------
86C
87C     INTLOG  - Logs messages
88C     JMALLOC - Dynamically allocate memory
89C     JFREE   - Free dynamically allocated memory
90C     JGETGG  - Reads the definition of a gaussian grid
91C     HGENGG  - Calculates original lat/long (before rotation) for
92C               a rotated gaussian grid
93C     HNEI12  - Finds neighbours for points for interpolation
94C     HWTS12  - Calculates weightings for points for interpolation
95C     FORCED_NEAREST_NEIGHBOUR - check forced interpolation method
96C
97C
98C     Reference
99C     ---------
100C
101C     None.
102C
103C
104C     Comments
105C     --------
106C
107C     None.
108C
109C
110C     Author
111C     ------
112C
113C     J.D.Chambers      ECMWF      February 2001
114C
115C
116C     Modifications
117C     -------------
118C
119C     None.
120C
121C----<
122C     -----------------------------------------------------------------|
123C*    Section 0. Definition of variables.
124C     -----------------------------------------------------------------|
125C
126      IMPLICIT NONE
127C
128#include "parim.h"
129#include "nifld.common"
130#include "nofld.common"
131C
132C     Parameters
133C
134      INTEGER JNORTH, JSOUTH, JWEST, JEAST, JW_E, JN_S, JLAT, JLON
135      INTEGER JP12PT, JP4PT, JPNEARN
136      PARAMETER (JP12PT  = 0)
137      PARAMETER (JP4PT   = 1)
138      PARAMETER (JPNEARN = 2)
139      PARAMETER (JNORTH = 1 )
140      PARAMETER (JWEST  = 2 )
141      PARAMETER (JSOUTH = 3 )
142      PARAMETER (JEAST  = 4 )
143      PARAMETER (JW_E  = 1 )
144      PARAMETER (JN_S  = 2 )
145      PARAMETER (JLAT  = 1 )
146      PARAMETER (JLON  = 2 )
147C
148C     Function arguments
149C
150      LOGICAL L12PNT
151      INTEGER KGNOLD, KGNNEW, KSIZE, NUMPTS
152      REAL AREA(4), POLE(2), OLDFLD(*), NEWFLD(KSIZE)
153      CHARACTER*1 HTYPE
154C
155C     Local variables
156C
157      INTEGER NEXT, LOOP, IRET, NLEN, NOPREV, NNPREV, NBYTES, NUMBER
158      INTEGER NEAREST
159      CHARACTER*1 NNTYPE
160C
161      LOGICAL LNEW, LFIRST, LVEGGY
162      INTEGER KSCHEME(1),NEIGH(12,1), KLA(1)
163      REAL PWTS(12,1)
164      POINTER (IPKSCHE, KSCHEME)
165      POINTER (IPNEIGH, NEIGH)
166      POINTER (IPKLA,   KLA)
167      POINTER (IPPWTS,  PWTS)
168C
169      REAL PDLO0(1),PDLO1(1),PDLO2(1),PDLO3(1),PDLAT(1)
170      POINTER (IPPDLO0, PDLO0)
171      POINTER (IPPDLO1, PDLO1)
172      POINTER (IPPDLO2, PDLO2)
173      POINTER (IPPDLO3, PDLO3)
174      POINTER (IPPDLAT, PDLAT)
175C
176      INTEGER IGG, IGGOLD, INN, INNOLD
177      INTEGER KPTS(1), KNPTS(1)
178      REAL GLATS(1), GNLATS(1)
179      INTEGER IOFFS(1)
180      POINTER (IPKPTS,  KPTS)
181      POINTER (IPKNPTS,  KNPTS)
182      POINTER (IPIOFFS, IOFFS)
183      POINTER (IPGLATS, GLATS)
184      POINTER (IPGNLATS, GNLATS)
185C
186      INTEGER ILL, ILLOLD
187      REAL RLAT(1),RLON(1)
188      POINTER (IPRLAT, RLAT)
189      POINTER (IPRLON, RLON)
190C
191      REAL OLD(1)
192      POINTER (IOLD,   OLD)
193C
194      DATA NOPREV/-1/, NNPREV/-1/, NNTYPE/' '/
195      DATA LNEW/.FALSE./, LFIRST/.TRUE./
196      DATA IGGOLD/-1/, INNOLD/-1/, ILLOLD/-1/, IOLD/-1/
197C
198      SAVE LNEW, LFIRST
199      SAVE IPKSCHE, IPNEIGH, IPKLA, IPPWTS
200      SAVE IPPDLO0, IPPDLO1, IPPDLO2, IPPDLO3, IPPDLAT
201      SAVE NOPREV,NNPREV,NNTYPE
202      SAVE IGGOLD,INNOLD,IPKPTS,IPKNPTS,IPIOFFS,IPGLATS,IPGNLATS
203      SAVE ILLOLD, IPRLAT, IPRLON, IOLD
204C
205C     Externals
206C
207      LOGICAL FORCED_NEAREST_NEIGHBOUR
208      INTEGER HNEI12, HGENGG
209#ifdef POINTER_64
210      INTEGER*8 JMALLOC
211#else
212      INTEGER JMALLOC
213#endif
214C
215C     -----------------------------------------------------------------|
216C     Section 1.  Initialise.
217C     -----------------------------------------------------------------|
218C
219  100 CONTINUE
220C
221      HRG2GG = 0
222C
223      CALL JDEBUG()
224
225      IF( L12PNT ) THEN
226        CALL INTLOG(JP_DEBUG,'HRG2GG: 12-pt interpolation',JPQUIET)
227      ELSE
228        CALL INTLOG(JP_DEBUG,'HRG2GG:  4-pt interpolation',JPQUIET)
229      ENDIF
230C
231      CALL CHKPREC()
232      IF( LPREC )THEN
233        CALL INTLOG(JP_DEBUG,
234     X   'HRG2GG: precipitation threshold applied',JPQUIET)
235      ELSE
236        CALL INTLOG(JP_DEBUG,
237     X   'HRG2GG: precipitation threshold not applied',JPQUIET)
238      ENDIF
239
240C     Use nearest neighbour if required
241      LVEGGY = FORCED_NEAREST_NEIGHBOUR(LMETHOD,NITABLE,NIPARAM)
242      IF( LVEGGY ) CALL INTLOG(JP_DEBUG,
243     X  'HRG2GG: nearest neighbour processing',JPQUIET)
244
245
246C
247C     Dynamically allocate memory for old gaussian grid information.
248C
249      IGG = KGNOLD*2
250C
251      IF( IGG.GT.IGGOLD ) THEN
252C
253        IF( IGGOLD.GT.0 ) CALL JFREE(IPKPTS)
254C
255        NBYTES = (IGG*JPRLEN) + (2*IGG+1)*JPBYTES
256C
257        IPKPTS = JMALLOC(NBYTES)
258#ifdef hpR64
259        IPKPTS = IPKPTS/(1024*1024*1024*4)
260#endif
261        IF( IPKPTS.EQ.0 ) THEN
262          CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET)
263          HRG2GG = 1
264          GOTO 900
265        ENDIF
266C
267        IPGLATS = IPKPTS  + (IGG*JPBYTES)
268        IPIOFFS = IPGLATS + (IGG*JPRLEN)
269C
270        IGGOLD = IGG
271        NOPREV = -1
272C
273      ENDIF
274C
275C     Build up offsets to start of each latitude in the original field.
276C
277      IF( KGNOLD.NE.NOPREV ) THEN
278        CALL JGETGG(KGNOLD,HTYPE,GLATS,KPTS,IRET)
279        IF( IRET.NE.0 ) THEN
280          CALL INTLOG(JP_ERROR,
281     X      'HRG2GG: JGETGG failed to get gaussian data',JPQUIET)
282          HRG2GG = 1
283          GOTO 900
284        ENDIF
285C
286        IOFFS(1) = 1
287        DO LOOP = 2, (KGNOLD*2+1)
288          IOFFS(LOOP) = IOFFS(LOOP-1) + KPTS(LOOP-1)
289        ENDDO
290C
291C       Allocate memory to hold the input field
292C       (in case OLDFLD and NEWFLD are the same arrays)
293C
294        IF( IOLD.GT.0 ) CALL JFREE(IOLD)
295C
296        NUMBER = (IOFFS(KGNOLD*2+1) - 1)
297        NBYTES = NUMBER * JPRLEN
298C
299        IOLD = JMALLOC(NBYTES)
300#ifdef hpR64
301        IOLD = IOLD/(1024*1024*1024*4)
302#endif
303        IF( IOLD.EQ.0 ) THEN
304          CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET)
305          HRG2GG = 1
306          GOTO 900
307        ENDIF
308C
309        NOPREV = KGNOLD
310      ENDIF
311C
312C     Preserve the input field
313C
314      NUMBER = (IOFFS(KGNOLD*2+1) - 1)
315      DO LOOP = 1, NUMBER
316        OLD(LOOP) = OLDFLD(LOOP)
317      ENDDO
318C
319C     -----------------------------------------------------------------|
320C     Section 2.  Generate the lat/long points for the output grid
321C     -----------------------------------------------------------------|
322C
323  200 CONTINUE
324C
325C
326C     Dynamically allocate memory for gaussian grid information.
327C
328      INN = KGNNEW*2
329C
330      IF( INN.GT.INNOLD ) THEN
331C
332        IF( INNOLD.GT.0 ) CALL JFREE(IPKNPTS)
333C
334        NBYTES = (INN*JPRLEN) + (INN*JPBYTES)
335C
336        IPKNPTS = JMALLOC(NBYTES)
337#ifdef hpR64
338        IPKNPTS = IPKNPTS/(1024*1024*1024*4)
339#endif
340        IF( IPKNPTS.EQ.0 ) THEN
341          CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET)
342          HRG2GG = 2
343          GOTO 900
344        ENDIF
345C
346        IPGNLATS = IPKNPTS  + (INN*JPBYTES)
347C
348        INNOLD = INN
349        NNPREV = -1
350      ENDIF
351C
352      IF( (KGNNEW.NE.NNPREV).OR.(HTYPE.NE.NNTYPE) ) THEN
353        IF( HTYPE.NE.'F' .AND. HTYPE.NE.'f' .AND.
354     X      HTYPE.NE.'O' .AND. HTYPE.NE.'o' ) THEN
355          HTYPE = 'R'
356        ENDIF
357        CALL JGETGG(KGNNEW,HTYPE,GNLATS,KNPTS,IRET)
358        IF( IRET.NE.0 ) THEN
359          CALL INTLOG(JP_ERROR,
360     X      'HRG2GG: JGETGG failed to get gaussian data',JPQUIET)
361          HRG2GG = 2
362          GOTO 900
363        ENDIF
364C
365        IF( HTYPE.EQ.'F' .OR. HTYPE.EQ.'f' ) THEN
366          NLEN = KGNNEW * KGNNEW * 8
367        ELSE
368          NLEN = 0
369          DO LOOP = 1, KGNNEW*2
370            NLEN = NLEN + KNPTS(LOOP)
371          ENDDO
372        ENDIF
373        NNPREV = KGNNEW
374        NNTYPE = HTYPE
375      ENDIF
376C
377C     Check that given array is big enough for the new field.
378C
379      IF( NLEN.GT.KSIZE ) THEN
380        CALL INTLOG(JP_ERROR,'HRG2GG: Given array size = ',KSIZE)
381        CALL INTLOG(JP_ERROR,'HRG2GG: Required size = ',NLEN)
382        HRG2GG = 2
383        GOTO 900
384      ENDIF
385C
386C     Dynamically allocate memory for lat/long arrays.
387C
388      ILL = NLEN
389      IF( ILL.GT.ILLOLD ) THEN
390C
391        LNEW = .TRUE.
392C
393        IF( ILLOLD.GT.0 ) CALL JFREE(IPRLON)
394C
395        NBYTES = 2*ILL*JPRLEN
396C
397        IPRLON = JMALLOC(NBYTES)
398#ifdef hpR64
399        IPRLON = IPRLON/(1024*1024*1024*4)
400#endif
401        IF( IPRLON.EQ.0 ) THEN
402          CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET)
403          HRG2GG = 2
404          GOTO 900
405        ENDIF
406C
407        IPRLAT = IPRLON + (ILL*JPRLEN)
408C
409        ILLOLD = ILL
410C
411      ENDIF
412C
413      IRET = HGENGG(AREA,POLE,KGNNEW,HTYPE,KNPTS,GNLATS,ILL,
414     X              RLAT,RLON,NUMPTS)
415      IF( IRET.NE.0 ) THEN
416        CALL INTLOG(JP_ERROR,
417     X    'HRG2GG: HGENGG failed to get lat/lon grid data',JPQUIET)
418        HRG2GG = 2
419        GOTO 900
420      ENDIF
421C
422C     -----------------------------------------------------------------|
423C     Section 3.  Find neighbours for each point for interpolation.
424C     -----------------------------------------------------------------|
425C
426  300 CONTINUE
427C
428C     Dynamically allocate memory for interpolation arrays.
429C
430      IF( LNEW ) THEN
431C
432        IF( .NOT.LFIRST ) CALL JFREE(IPPDLO0)
433C
434        NBYTES = (17*JPRLEN + 14*JPBYTES) * ILL
435C
436        IPPDLO0 = JMALLOC(NBYTES)
437#ifdef hpR64
438        IPPDLO0 = IPPDLO0/(1024*1024*1024*4)
439#endif
440        IF( IPPDLO0.EQ.0 ) THEN
441          CALL INTLOG(JP_ERROR,'HRG2GG: Memory allocation fail',JPQUIET)
442          HRG2GG = 3
443          GOTO 900
444        ENDIF
445C
446        IPPDLO1 = IPPDLO0 + (ILL*JPRLEN)
447        IPPDLO2 = IPPDLO1 + (ILL*JPRLEN)
448        IPPDLO3 = IPPDLO2 + (ILL*JPRLEN)
449        IPPDLAT = IPPDLO3 + (ILL*JPRLEN)
450        IPPWTS  = IPPDLAT + (ILL*JPRLEN)
451        IPKSCHE = IPPWTS  + (12*ILL*JPRLEN)
452        IPKLA   = IPKSCHE + (ILL*JPBYTES)
453        IPNEIGH = IPKLA   + (ILL*JPBYTES)
454C
455        LFIRST = .FALSE.
456        LNEW   = .FALSE.
457C
458      ENDIF
459C
460C     Find neighbours.
461C
462      IRET = HNEI12(L12PNT,NLEN,RLAT,RLON,KGNOLD,KPTS,GLATS,
463     X              KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH)
464      IF( IRET.NE.0 ) THEN
465        CALL INTLOG(JP_ERROR,
466     X    'HRG2GG: HNEI12 failed to find neighbours',JPQUIET)
467        HRG2GG = 3
468        GOTO 900
469      ENDIF
470C
471C     -----------------------------------------------------------------|
472C     Section 4.  Perform the 12-point horizontal interpolation.
473C     -----------------------------------------------------------------|
474C
475  400 CONTINUE
476C
477C     Setup the 12-point horizontal interpolation weights
478C
479      CALL HWTS12
480     X  (NLEN,KSCHEME,KLA,PDLAT,GLATS,PDLO0,PDLO1,PDLO2,PDLO3,NEIGH,
481     X   PWTS)
482C
483C     Calculate the interpolated grid point values
484C
485      DO LOOP = 1, NLEN
486       IF( LVEGGY) THEN
487            NEAREST = 1
488            IF( PWTS( 2,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 2
489            IF( PWTS( 3,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 3
490            IF( PWTS( 4,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 4
491            IF( PWTS( 5,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 5
492            IF( PWTS( 6,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 6
493            IF( PWTS( 7,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 7
494            IF( PWTS( 8,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 8
495            IF( PWTS( 9,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST = 9
496            IF( PWTS(10,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =10
497            IF( PWTS(11,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =11
498            IF( PWTS(12,LOOP).GT.PWTS(NEAREST,LOOP) ) NEAREST =12
499            NEWFLD(LOOP) = OLD(NEIGH( NEAREST,LOOP))
500        ELSE
501          IF( KSCHEME(LOOP).EQ.JP12PT ) THEN
502            NEWFLD(LOOP) =
503     X        OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) +
504     X        OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) +
505     X        OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) +
506     X        OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP) +
507     X        OLD(NEIGH( 5,LOOP)) * PWTS( 5,LOOP) +
508     X        OLD(NEIGH( 6,LOOP)) * PWTS( 6,LOOP) +
509     X        OLD(NEIGH( 7,LOOP)) * PWTS( 7,LOOP) +
510     X        OLD(NEIGH( 8,LOOP)) * PWTS( 8,LOOP) +
511     X        OLD(NEIGH( 9,LOOP)) * PWTS( 9,LOOP) +
512     X        OLD(NEIGH(10,LOOP)) * PWTS(10,LOOP) +
513     X        OLD(NEIGH(11,LOOP)) * PWTS(11,LOOP) +
514     X        OLD(NEIGH(12,LOOP)) * PWTS(12,LOOP)
515C
516          ELSE IF( KSCHEME(LOOP).EQ.JP4PT ) THEN
517            NEWFLD(LOOP) =
518     X        OLD(NEIGH( 1,LOOP)) * PWTS( 1,LOOP) +
519     X        OLD(NEIGH( 2,LOOP)) * PWTS( 2,LOOP) +
520     X        OLD(NEIGH( 3,LOOP)) * PWTS( 3,LOOP) +
521     X        OLD(NEIGH( 4,LOOP)) * PWTS( 4,LOOP)
522C
523          ELSE
524            DO NEXT = 1, 4
525              IF( NEIGH(NEXT,LOOP).NE.0 )
526     X          NEWFLD(LOOP) = OLD(NEIGH(NEXT,LOOP))
527            ENDDO
528C
529         ENDIF
530        ENDIF
531      ENDDO
532C
533C     -----------------------------------------------------------------|
534C     Section 9.  Return.
535C     -----------------------------------------------------------------|
536C
537  900 CONTINUE
538C
539      RETURN
540      END
541
542