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