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 HNEI12
12     X  (L12PNT,KLEN,RLAT,RLON,KGAUSS,KPTS,PLATIN,
13     X   KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH)
14C
15C---->
16C**** HNEI12
17C
18C     Purpose
19C     -------
20C
21C     This routine accepts a vector of points and finds the 12
22C     neighbours for each point suitable for horizontal interpolation.
23C
24C
25C     Interface
26C     ---------
27C
28C     IRET = HNEI12(L12PNT,KLEN,RLAT,RLON,KGAUSS,KPTS,PLATIN,
29C    X              KSCHEME,PDLAT,PDLO0,PDLO1,PDLO2,PDLO3,KLA,NEIGH)
30C
31C
32C     Input parameters
33C     ----------------
34C
35C     L12PNT - Chooses between 12-point and 4-point interpolation
36C              = .TRUE. for 12-point horizontal
37C              = .FALSE. for 4-point
38C     KLEN   - Number of points along the vector
39C     RLAT   - List of latitudes for points.
40C     RLON   - List of longitudes for points.
41C     KGAUSS - Original gaussian field number.
42C     KPTS   - List giving number of points in each latitude in original
43C              gaussian field.
44C     PLATIN - Original gaussian field latitudes (north and south).
45C
46C
47C     Output parameters
48C     -----------------
49C
50C     KSCHEME - Flag showing interpolation scheme to use for point
51C               0 = 12-point
52C               1 = 4-point bilinear
53C               2 = nearest neighbour
54C     PDLAT   - Meridian linear weight.
55C     PDLO0   - Zonal linear weight for the latitude of point 5.
56C     PDLO1   - Zonal linear weight for the latitude of point 1.
57C     PDLO2   - Zonal linear weight for the latitude of point 3.
58C     PDLO3   - Zonal linear weight for the latitude of point 11.
59C     KLA     - Latitude number in original field of latitude north of
60C               each point in the vector.
61C     NEIGH   - List of indices in the original field of neighbouring
62C               point values for each point in the vector.
63C
64C     Returns 0 if function successful, non-zero otherwise.
65C
66C     Common block usage
67C     ------------------
68C
69C     None.
70C
71C
72C     Method
73C     ------
74C
75C     Numbering of the points (I is the interpolation point):
76C
77C                   13       5       6      14
78C
79C                    7       1       2       8
80C                               (I)
81C                    9       3       4      10
82C
83C                   15      11      12      16
84C
85C     The 12-point interpolation is not possible if either of the top
86C     two rows is above the original field northern latitude. The
87C     nearest neighbour is used if both rows are above, and a 4-pt
88C     bilinear interpolation is used if the top row is above.
89C     Similarily, if either of the bottom two rows is below the original
90C     field southern latitude.
91C
92C
93C     Externals
94C     ---------
95C
96C     INTLOG  - Log error message.
97C     JNORSGG - Finds gaussian latitude north or south of given latitude
98C     JDEBUG  - Tests if debug output required.
99C     JMALLOC - Dynamically allocate memory
100C     JFREE   - Free dynamically allocated memory
101C
102C
103C     Reference
104C     ---------
105C
106C     ECMWF Meteorological Bulletin M1.6/7
107C     IFS Documentation
108C     Part VI: Technical and Computational Procedures (CY21R4)
109C     March 2000
110C     Section 2.3
111C
112C
113C     Comments
114C     --------
115C
116C     None.
117C
118C
119C     Author
120C     ------
121C
122C     J.D.Chambers      ECMWF      January 2001
123C
124C
125C     Modifications
126C     -------------
127C
128C     None.
129C
130C----<
131C     -----------------------------------------------------------------|
132C*    Section 0. Definition of variables.
133C     -----------------------------------------------------------------|
134C
135      IMPLICIT NONE
136#include "parim.h"
137C
138C     Parameters
139C
140      INTEGER JP12PT, JP4PT, JPNEARN
141      PARAMETER (JP12PT  = 0)
142      PARAMETER (JP4PT   = 1)
143      PARAMETER (JPNEARN = 2)
144C
145C     Function arguments
146C
147      LOGICAL L12PNT
148      INTEGER KLEN, KGAUSS, KPTS(KGAUSS*2)
149      REAL RLAT(KLEN), RLON(KLEN), PLATIN(KGAUSS*2)
150C
151      INTEGER KLA(KLEN), KSCHEME(KLEN)
152      REAL PDLAT(KLEN),PDLO0(KLEN),PDLO1(KLEN),PDLO2(KLEN),PDLO3(KLEN)
153      INTEGER NEIGH(12,*)
154C
155C     Local variables
156C
157      INTEGER NEXT, LOOP, NBYTES
158      INTEGER NLAT0, NLAT1, NLAT2, NLAT3, NLAT(12), NLONG(12)
159      REAL STEP0, STEP1, STEP2, STEP3
160      INTEGER NPREV, ISIZE, ISIZOLD, IOFFSET(1)
161      POINTER( IPOFF, IOFFSET )
162C
163      DATA ISIZOLD/-1/, NPREV/-1/
164      SAVE ISIZOLD,IPOFF,NPREV
165C
166C     Externals
167C
168      INTEGER JNORSGG
169#ifdef POINTER_64
170      INTEGER*8 JMALLOC
171#else
172      INTEGER JMALLOC
173#endif
174C
175C     -----------------------------------------------------------------|
176C     Section 1.  Initialise.
177C     -----------------------------------------------------------------|
178C
179  100 CONTINUE
180C
181      HNEI12 = 0
182C
183      CALL JDEBUG()
184C
185C     Get memory for offset array.
186C
187      ISIZE = KGAUSS*2 + 1
188      IF( ISIZE.GT.ISIZOLD ) THEN
189C
190        IF( ISIZOLD.GT.0 ) CALL JFREE(IPOFF)
191C
192        NBYTES = ISIZE * JPBYTES
193        IPOFF = JMALLOC(NBYTES)
194#ifdef hpR64
195        IPOFF = IPOFF/(1024*1024*1024*4)
196#endif
197        IF( IPOFF.EQ.0 ) THEN
198          CALL INTLOG(JP_WARN,'HNEI12: Memory allocate fail ',JPQUIET)
199          HNEI12 = 1
200          GOTO 900
201        ENDIF
202C
203        ISIZOLD = ISIZE
204        NPREV = -1
205C
206      ENDIF
207C
208C     Build up offsets to start of each latitude in the original field.
209C
210      IF( KGAUSS.NE.NPREV ) THEN
211        IOFFSET(1) = 1
212        DO LOOP = 2, (KGAUSS*2+1)
213          IOFFSET(LOOP) = IOFFSET(LOOP-1) + KPTS(LOOP-1)
214        ENDDO
215        NPREV = KGAUSS
216      ENDIF
217C
218C     -----------------------------------------------------------------|
219C       Section 2.  Loop through points along the vector.
220C     -----------------------------------------------------------------|
221C
222  200 CONTINUE
223C
224      DO NEXT = 1,KLEN
225C
226C       Clear current neighbour indices
227C
228        DO LOOP = 1, 12
229          NEIGH(LOOP,NEXT) = 0
230        ENDDO
231C
232C       Find latitude numbers for neighbours.
233C
234C       First find latitude numbers to north and south of the point
235C
236C     -----------------------------------------------------------------|
237C       Section 3. Nearest neighbour to be used.
238C     -----------------------------------------------------------------|
239C
240  300   CONTINUE
241C
242        IF( RLAT(NEXT).GT.PLATIN(1) ) THEN
243          CALL INTLOG(JP_DEBUG,
244     X      'HNEI12: Nearest neighbour to south used for point ',NEXT)
245C
246C         Above northern latitude of original field, so use nearest
247C         neighbour (point 3 or 4)
248C
249          STEP2 = 360.0 / REAL(KPTS(1))
250          NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2)
251          IF( NLONG( 3).GT.KPTS(1) ) NLONG( 3) = 1
252          NLONG( 4) = NLONG(3) + 1
253          IF( NLONG( 4).GT.KPTS(1) ) NLONG( 4) = 1
254C
255          IF( ABS((RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 )) .LE.
256     X        ABS((RLON(NEXT) - REAL(NLONG(4)-1)*STEP2 )) ) THEN
257            NEIGH(3,NEXT) = NLONG(3)
258          ELSE
259            NEIGH(4,NEXT) = NLONG(4)
260          ENDIF
261          KSCHEME(NEXT) = JPNEARN
262          GOTO 600
263        ENDIF
264C
265        IF( RLAT(NEXT).LT.PLATIN(KGAUSS*2) ) THEN
266          CALL INTLOG(JP_DEBUG,
267     X      'HNEI12: Nearest neighbour to north used for point ',NEXT)
268C
269C         Below southern latitude of original field, so use nearest
270C         neighbour (point 1 or 2)
271C
272          STEP1 = 360.0 / REAL(KPTS(KGAUSS*2))
273          NLONG( 1) = 1 + INT(RLON(NEXT)/STEP1)
274          IF( NLONG( 1).GT.KPTS(KGAUSS*2) ) NLONG( 1) = 1
275          NLONG( 2) = NLONG(1) + 1
276          IF( NLONG( 2).GT.KPTS(KGAUSS*2) ) NLONG( 2) = 1
277C
278          IF( ABS((RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 )) .LE.
279     X        ABS((RLON(NEXT) - REAL(NLONG(2)-1)*STEP1 )) ) THEN
280            NEIGH(1,NEXT) = IOFFSET(KGAUSS*2) + NLONG(1) - 1
281          ELSE
282            NEIGH(2,NEXT) = IOFFSET(KGAUSS*2) + NLONG(2) - 1
283          ENDIF
284          KSCHEME(NEXT) = JPNEARN
285          GOTO 600
286        ENDIF
287C
288        NLAT1 = JNORSGG(RLAT(NEXT),PLATIN,KGAUSS,1)
289        NLAT2 = JNORSGG(RLAT(NEXT),PLATIN,KGAUSS,0)
290C
291        IF( NLAT2.EQ.NLAT1) NLAT2 = NLAT1 + 1
292C
293C       Find outer latitude numbers to north and south of the point
294C
295        NLAT0 = NLAT1 - 1
296        NLAT3 = NLAT2 + 1
297C
298C     -----------------------------------------------------------------|
299C       Section 4. 4-point bilinear to be used.
300C     -----------------------------------------------------------------|
301C
302  400   CONTINUE
303C
304        IF( (NLAT0.EQ.0).OR.(NLAT3.GT.KGAUSS*2).OR.(.NOT.L12PNT) ) THEN
305Cjdc      CALL INTLOG(JP_DEBUG,
306Cjdc X      'HNEI12: 4-pt bilinear interpolation used for point ',NEXT)
307C
308C         Between two northern latitudes or two southern latitudes of
309C         original field, so use 4-point bilinear interpolation
310C         (points 1, 2, 3 and 4)
311C
312C         Points 1 and 2
313C
314          STEP1 = 360.0 / REAL(KPTS(NLAT1))
315          NLONG(1) = 1 + INT(RLON(NEXT)/STEP1)
316          IF( NLONG(1).GT.KPTS(NLAT1) ) NLONG(1) = 1
317          NLONG(2) = NLONG(1) + 1
318          IF( NLONG(2).GT.KPTS(NLAT1) ) NLONG(2) = 1
319C
320          NEIGH(1,NEXT) = IOFFSET(NLAT1) + NLONG(1) - 1
321          NEIGH(2,NEXT) = IOFFSET(NLAT1) + NLONG(2) - 1
322C
323C         Points 3 and 4
324C
325          STEP2 = 360.0 / REAL(KPTS(NLAT2))
326          NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2)
327          IF( NLONG( 3).GT.KPTS(NLAT2) ) NLONG( 3) = 1
328          NLONG( 4) = NLONG(3) + 1
329          IF( NLONG( 4).GT.KPTS(NLAT2) ) NLONG( 4) = 1
330C
331          NEIGH(3,NEXT) = IOFFSET(NLAT2) + NLONG(3) - 1
332          NEIGH(4,NEXT) = IOFFSET(NLAT2) + NLONG(4) - 1
333C
334C         Calculate zonal linear weights for the latitudes.
335C
336          PDLO1(NEXT) = ( RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 ) / STEP1
337          PDLO2(NEXT) = ( RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 ) / STEP2
338C
339          KLA(NEXT) = NLAT1
340          PDLAT(NEXT) = ( RLAT(NEXT) - PLATIN(NLAT1) ) /
341     X                  ( PLATIN(NLAT2) - PLATIN(NLAT1) )
342C
343          KSCHEME(NEXT) = JP4PT
344          GOTO 600
345        ENDIF
346C
347C     -----------------------------------------------------------------|
348C       Section 5. 12-point interpolation to be used.
349C     -----------------------------------------------------------------|
350C
351  500   CONTINUE
352C
353        KSCHEME(NEXT) = JP12PT
354C
355        KLA(NEXT) = NLAT1
356        PDLAT(NEXT) = ( RLAT(NEXT) - PLATIN(NLAT1) ) /
357     X                ( PLATIN(NLAT2) - PLATIN(NLAT1) )
358C
359C       Northernmost
360C
361        NLAT(5) = NLAT0
362        NLAT(6) = NLAT0
363C
364C       Northern
365C
366        NLAT(7) = NLAT1
367        NLAT(1) = NLAT1
368        NLAT(2) = NLAT1
369        NLAT(8) = NLAT1
370C
371C       Southern
372C
373        NLAT( 9) = NLAT2
374        NLAT( 3) = NLAT2
375        NLAT( 4) = NLAT2
376        NLAT(10) = NLAT2
377C
378C       Southernmost
379C
380        NLAT(11) = NLAT3
381        NLAT(12) = NLAT3
382C
383C     -----------------------------------------------------------------|
384C       Find longitude numbers for neighbours.
385C     -----------------------------------------------------------------|
386C
387C       Northernmost
388C
389        STEP0 = 360.0 / REAL(KPTS(NLAT0))
390        NLONG(5) = 1 + INT(RLON(NEXT)/STEP0)
391        IF( NLONG(5).GT.KPTS(NLAT0) ) NLONG(5) = 1
392        NLONG(6) = NLONG(5) + 1
393        IF( NLONG(6).GT.KPTS(NLAT0) ) NLONG(6) = 1
394C
395C       Northern
396C
397        STEP1 = 360.0 / REAL(KPTS(NLAT1))
398        NLONG(1) = 1 + INT(RLON(NEXT)/STEP1)
399        IF( NLONG(1).GT.KPTS(NLAT1) ) NLONG(1) = 1
400        NLONG(7) = NLONG(1) - 1
401        IF( NLONG(7).LT.1 ) NLONG(7) = KPTS(NLAT1)
402        NLONG(2) = NLONG(1) + 1
403        IF( NLONG(2).GT.KPTS(NLAT1) ) NLONG(2) = 1
404        NLONG(8) = NLONG(2) + 1
405        IF( NLONG(8).GT.KPTS(NLAT1) ) NLONG(8) = 1
406C
407C       Southern
408C
409        STEP2 = 360.0 / REAL(KPTS(NLAT2))
410        NLONG( 3) = 1 + INT(RLON(NEXT)/STEP2)
411        IF( NLONG( 3).GT.KPTS(NLAT2) ) NLONG( 3) = 1
412        NLONG( 9) = NLONG(3) - 1
413        IF( NLONG( 9).LT.1 ) NLONG( 9) = KPTS(NLAT2)
414        NLONG( 4) = NLONG(3) + 1
415        IF( NLONG( 4).GT.KPTS(NLAT2) ) NLONG( 4) = 1
416        NLONG(10) = NLONG(4) + 1
417        IF( NLONG(10).GT.KPTS(NLAT2) ) NLONG(10) = 1
418C
419C       Southernmost
420C
421        STEP3 = 360.0 / REAL(KPTS(NLAT3))
422        NLONG(11) = 1 + INT(RLON(NEXT)/STEP3)
423        IF( NLONG(11).GT.KPTS(NLAT3) ) NLONG(11) = 1
424        NLONG(12) = NLONG(11) + 1
425        IF( NLONG(12).GT.KPTS(NLAT3) ) NLONG(12) = 1
426C
427C     -----------------------------------------------------------------|
428C       Calculate zonal linear weights for the latitudes.
429C     -----------------------------------------------------------------|
430C
431        PDLO0(NEXT) = ( RLON(NEXT) - REAL(NLONG(5)-1)*STEP0 ) / STEP0
432C
433        PDLO1(NEXT) = ( RLON(NEXT) - REAL(NLONG(1)-1)*STEP1 ) / STEP1
434C
435        PDLO2(NEXT) = ( RLON(NEXT) - REAL(NLONG(3)-1)*STEP2 ) / STEP2
436C
437        PDLO3(NEXT) = ( RLON(NEXT) - REAL(NLONG(11)-1)*STEP3 ) / STEP3
438C
439C     -----------------------------------------------------------------|
440C       Store indices of the neighbours.
441C     -----------------------------------------------------------------|
442C
443        DO LOOP = 1, 12
444          NEIGH(LOOP,NEXT) = IOFFSET(NLAT(LOOP)) + NLONG(LOOP) - 1
445        ENDDO
446C
447C     -----------------------------------------------------------------|
448C     Section 6.  End of loop along vector of points.
449C     -----------------------------------------------------------------|
450C
451  600   CONTINUE
452C
453      ENDDO
454C
455C     -----------------------------------------------------------------|
456C     Section 9.  Return
457C     -----------------------------------------------------------------|
458C
459  900 CONTINUE
460C
461      RETURN
462      END
463