1C
2C Copyright 1981-2016 ECMWF.
3C
4C This software is licensed under the terms of the Apache Licence
5C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6C
7C In applying this licence, ECMWF does not waive the privileges and immunities
8C granted to it by virtue of its status as an intergovernmental organisation
9C nor does it submit to any jurisdiction.
10C
11
12      INTEGER FUNCTION WAVEXX2(NPARAM, NUMLATS, NPTS,
13     X                         NLATS, STEPNS, STEPWE,
14     X                         OLDWAVE, NEWWAVE, NORTH, WEST, PMISS)
15C
16C---->
17C*****WAVEXX2*
18C
19C     PURPOSE
20C     -------
21C
22C     Interpolates wave fields (except 2D spectra).
23C
24C
25C     INTERFACE
26C     ---------
27C
28C     IRET = WAVEXX2(NPARAM, NUMLATS, NPTS,
29C    X               NLATS, STEPNS, STEPWE,
30C    X               OLDWAVE, NEWWAVE, NORTH, WEST, PMISS)
31C
32C     Input arguments
33C     ---------------
34C
35C     NPARAM  - Parameter ID
36C     NUMLATS - Input lat number north-south
37C     NPTS    - Array giving number of points along each latitude
38C               for input field (empty latitudes have entry 0)
39C     NLATS   - Number of points N-S in new grid
40C     STEPNS  - Output grid north-south resolution (degrees)
41C     STEPWE  - Output grid west-east resolution (degrees)
42C     OLDWAVE - Original wave field
43C     NORTH   - Output grid northernmost latitude (degrees)
44C     WEST    - Output grid westernmost longitude (degrees)
45C     PMISS   - Missing data value
46C
47C     Output arguments
48C     ----------------
49C
50C     NEWWAVE - New wave field
51C
52C     Function returns 0 if the interpolation was OK.
53C
54C
55C     METHOD
56C     ------
57C
58C     Builds the index of neighbouring points for the output grid.
59C     Then works through the output grid points, checking for subarea
60C     boundaries and looking up neighbouring point values and weights
61C     (which may be missing data).
62C
63C
64C     EXTERNALS
65C     ---------
66C
67C     WAVEIDX - Determines which nearest neighbour points to use for
68C               interpolating to new output grid point
69C     NUMPTWE - Calculates number of grid points between west/east
70C               area boundaries
71C     INTLOG  - Log error message
72C     JMALLOC - Dynamically allocate memory (deallocation -JFREE- not used)
73C
74C
75C     REFERENCE
76C     ---------
77C
78C     None.
79C
80C
81C     Author.
82C     -------
83C
84C     S. Curic      ECMWF    Jun 2009
85C
86C
87C
88C----<
89C
90      IMPLICIT NONE
91C
92C     Parameters
93C
94#include "parim.h"
95C
96      INTEGER JPROUTINE, JPMXLAT,
97     X        JPNW, JPNE, JPSW, JPSE, JPN, JPS,
98     X        JPDISNW, JPDISNE, JPDISSW, JPDISSE
99      PARAMETER( JPROUTINE = 40100 )
100      PARAMETER( JPMXLAT   =  1801 ) ! allow up to 0.1 degree resolution
101      PARAMETER( JPNW=1, JPNE=2, JPSW=3, JPSE=4, JPN=5, JPS=6,
102     X           JPDISNW=7, JPDISNE=8, JPDISSW=9, JPDISSE=10 )
103C
104#include "nifld.common"
105#include "nofld.common"
106#include "grspace.h"
107C
108C     Arguments
109C
110      INTEGER NPARAM
111      INTEGER NUMLATS, NPTS(*), NLATS
112      REAL STEPNS, STEPWE, OLDWAVE(*), NEWWAVE(*), NORTH, WEST, PMISS
113C
114C     Local variables
115C
116      INTEGER IEOFSET, INDEX, ISTART, IWEST, IWOFSET, KNEWNUM,
117     X        KOLDNUM, LOOP, MISSLAT, NCOL, NEXT,
118     X        NEXTWV, NROW, NUMNEW(JPMXLAT)
119      REAL    AWEST, EAST, NEWLATS(JPMXLAT), OLDLATS(JPMXLAT),
120     X        ONORTH, OSOUTH, PTLAT, PTLONG, RLATINC, ROWINC,
121     X        SOUTH, OEAST, OWEST
122      INTEGER INE, INW, ISE, ISW
123      LOGICAL LDIREC
124      REAL    C1, C2, CC, CNE_PT, CNW_PT, CSE_PT, CSW_PT, DI1N, DI1S,
125     X        DI2N, DI2S, DISNE, DISNW, DISSE, DISSW, DK1, DK2, NE_PT,
126     X        NW_PT, RAD, S1, S2, SE_PT, SNE_PT, SNW_PT, SS, SSE_PT,
127     X        SSW_PT, SW_PT, U1, U2
128C
129      DATA RAD/0.01745329251/
130C
131C     Large, resolution-dependent arrays
132C
133      LOGICAL IS_ALLOCATED_LARGEARRAY
134      DATA    IS_ALLOCATED_LARGEARRAY /.FALSE./
135      SAVE    IS_ALLOCATED_LARGEARRAY
136C
137      INTEGER NEWIDX(4,JPARRAYDIM_WAVE)
138#ifndef _CRAYFTN
139#ifdef POINTER_64
140      INTEGER*8 P_NEWIDX
141#endif
142#endif
143      POINTER(P_NEWIDX,NEWIDX)
144      SAVE P_NEWIDX
145C
146      REAL*4 DISTNEW(10,JPARRAYDIM_WAVE)
147#ifndef _CRAYFTN
148#ifdef POINTER_64
149      INTEGER*8 P_DISTNEW
150#endif
151#endif
152      POINTER(P_DISTNEW,DISTNEW)
153      SAVE P_DISTNEW
154C
155C     Externals
156C
157      INTEGER  WAVEIDX, NUMPTWE
158#ifdef POINTER_64
159      INTEGER*8 JMALLOC
160#else
161      INTEGER JMALLOC
162#endif
163      LOGICAL IS_WAVE_DIRECTION
164      EXTERNAL WAVEIDX, NUMPTWE, JMALLOC, IS_WAVE_DIRECTION
165C
166C     -----------------------------------------------------------------|
167C*    Section 1. Initalisation.
168C     -----------------------------------------------------------------|
169C
170  100 CONTINUE
171C
172      WAVEXX2 = 0
173C
174C     Initialise the bitmap value lookup function
175C
176      MISSLAT = 0
177
178      ONORTH = FLOAT(NIAREA(1))/PPMULT
179      OSOUTH = FLOAT(NIAREA(3))/PPMULT
180      RLATINC = FLOAT(NIGRID(2))/PPMULT
181C
182C     Calculate number of latitudes if grid had been full from
183C     North pole to South pole
184C
185      IF( NUMLATS.GT.JPMXLAT ) THEN
186        CALL INTLOG(JP_ERROR,
187     X    'WAVEXX2: Number of latitudes in input grid = ',NUMLATS)
188        CALL INTLOG(JP_ERROR,
189     X    'WAVEXX2: And is greater than allowed maximum = ',JPMXLAT)
190        WAVEXX2 = JPROUTINE + 1
191        GOTO 900
192      ENDIF
193C
194C
195C     Fill an array with the number of points at each latitude for the
196C     input field.
197C
198      IF(NIREPR.EQ.JPREDLL) THEN
199C
200C       Input field is a reduced latitude/longitude grid
201C
202C       .. but it may be 'pseudo-gaussian' in layout
203C       (ie global, symmetric about the equator but no latitude
204C        at the equator)
205C
206        IF( (ONORTH.NE.90.0).AND.(MOD(NUMLATS,2).EQ.0) ) THEN
207C
208          ONORTH = FLOAT(NIAREA(1))
209          OSOUTH = FLOAT(NIAREA(3))
210          RLATINC = FLOAT(NIGRID(2))
211
212        ENDIF
213C
214        DO LOOP = 1, NUMLATS
215          OLDLATS(LOOP) = 90.0 - (LOOP-1)*RLATINC
216        ENDDO
217C
218      ELSE
219C
220C       Input field is a regular latitude/longitude grid
221C
222        DO LOOP = 1, NUMLATS
223          OLDLATS(LOOP) = ONORTH - (LOOP-1)*RLATINC
224        ENDDO
225C
226        NPTS(1:NUMLATS) = NIWE
227C
228      ENDIF
229C
230C     Allocate large, resolution-dependent arrays
231      IF (.NOT. IS_ALLOCATED_LARGEARRAY) THEN
232        P_NEWIDX  = JMALLOC( 4*JPARRAYDIM_WAVE*JPBYTES)  !  4*len*bytes
233        P_DISTNEW = JMALLOC(10*JPARRAYDIM_WAVE*4)        ! 10*len*bytes
234        IF( (P_NEWIDX.EQ.0) .OR. (P_DISTNEW.EQ.0) ) THEN
235          CALL INTLOG(JP_ERROR,'WAVEXX2: JMALLOC fail',JPQUIET)
236          WAVEXX2 = JPROUTINE + 1
237          GOTO 900
238        ENDIF
239        IS_ALLOCATED_LARGEARRAY = .TRUE.
240      END IF
241C
242C     -----------------------------------------------------------------|
243C*    Section 2. Setup number of points at each latitude for the
244C                output latitude/longitude field.
245C     -----------------------------------------------------------------|
246C
247  200 CONTINUE
248C
249      IF( (NOREPR.EQ.JPQUASI) .OR. (NOREPR.EQ.JPGAUSSIAN) ) THEN
250C
251C       Reduced (quasi-regular) gaussian output
252C
253        KNEWNUM = NOGAUSS*2
254        DO LOOP = 1, KNEWNUM
255          NUMNEW(LOOP)  = NOLPTS(LOOP)
256          NEWLATS(LOOP) = ROGAUSS(LOOP)
257        ENDDO
258C
259      ELSE IF( NOREPR.EQ.JPREDLL ) THEN
260C
261C       Reduced (quasi-regular) lat/long output
262C
263        KNEWNUM = NOREDLL
264        DO LOOP = 1, KNEWNUM
265          NUMNEW(LOOP)  = NOLPTS(LOOP)
266          NEWLATS(LOOP) = ROREDLL(LOOP)
267        ENDDO
268C
269      ELSE
270C
271C       Regular output
272C
273        MISSLAT = NINT((90.0 - NORTH)/STEPNS)
274        DO LOOP = 1, MISSLAT
275          NUMNEW(LOOP)    = 0
276        ENDDO
277        DO LOOP = 1, NLATS
278          NUMNEW(LOOP+MISSLAT) = NINT(360.0/STEPWE)
279        ENDDO
280C
281        KNEWNUM = MISSLAT + NLATS
282        DO LOOP = 1, KNEWNUM
283          NEWLATS(LOOP) = 90.0 - (LOOP-1)*STEPNS
284        ENDDO
285C
286      ENDIF
287C
288C     -----------------------------------------------------------------|
289C*    Section 3. Get the input GRIB field characteristics.
290C     -----------------------------------------------------------------|
291C
292  300 CONTINUE
293C
294C     Calculate the indices of the input grid points to be used for
295C     the output points
296C
297      OWEST = FLOAT(NIAREA(2))/PPMULT
298      OEAST = FLOAT(NIAREA(4))/PPMULT
299
300      WAVEXX2 = WAVEIDX(NUMLATS,NPTS,OLDLATS,OWEST,OEAST,
301     X                  KNEWNUM, NUMNEW, NEWLATS, NEWIDX, DISTNEW)
302      IF( WAVEXX2.NE.0 ) THEN
303        CALL INTLOG(JP_ERROR,
304     X    'WAVEXX2: Unable to calculate output grid indices',JPQUIET)
305        WAVEXX2 = JPROUTINE + 3
306        GOTO 900
307      ENDIF
308C
309C     Wave direction parameters need special handling
310      LDIREC = IS_WAVE_DIRECTION(NPARAM)
311C
312C     -----------------------------------------------------------------|
313C*    Section 4. Work through the output subarea.
314C     -----------------------------------------------------------------|
315C
316  400 CONTINUE
317C
318C     Fill in the wave spectra values
319C
320      NEXT = 0
321      NEXTWV = 0
322C
323      SOUTH = NOAREA(3)/PPMULT
324      EAST  = NOAREA(4)/PPMULT
325      ISTART = 0
326C
327C     Work down through latitudes from north to south.
328C
329      DO NROW = 1, KNEWNUM
330C
331C       If inside north-south (subarea) boundaries ..
332C
333        IF( (NOREPR.EQ.JPGAUSSIAN).OR.(NOREPR.EQ.JPQUASI) ) THEN
334          PTLAT = ROGAUSS(NROW)
335        ELSE
336          PTLAT = 90.0 - (NROW-1)*STEPNS
337        ENDIF
338C
339        IF( (PTLAT.LE.NORTH).AND.(ABS(PTLAT-SOUTH).GT.-0.0005) ) THEN
340C
341C         Calculate number of points between west boundary of area and
342C         Greenwich
343C
344          ROWINC = 360.0/NUMNEW(NROW)
345C
346          IWEST = INT(WEST/ROWINC)
347          IF (ABS(WEST - IWEST*ROWINC).GT.0.000001) THEN
348            IF (WEST.GT.0.0) IWEST = IWEST + 1
349            IF (WEST.LT.0.0) IWEST = IWEST - 1
350          ENDIF
351          AWEST = IWEST * ROWINC
352          IWOFSET = NUMPTWE(AWEST,0.0,ROWINC)
353          IEOFSET = NUMPTWE(AWEST,EAST,ROWINC)
354C
355C         Work through subarea longitudes from west to east.
356C
357          DO NCOL = 1, NUMNEW(NROW)
358            PTLONG = AWEST + (NCOL-1)*ROWINC
359            NEXT = NUMPTWE(AWEST,PTLONG,ROWINC)
360            IF( (NEXT.LE.IEOFSET).AND.(NEXT.GE.0) ) THEN
361C
362C             .. and inside west-east (subarea) boundaries
363C
364              NEXT = 1 + NEXT - IWOFSET
365              IF( NEXT.LE.0) NEXT = NEXT + NUMNEW(NROW)
366              NEXT = NEXT + ISTART
367              NEXTWV = NEXTWV + 1
368C
369              INW = NEWIDX(JPNW,NEXT)
370              INE = NEWIDX(JPNE,NEXT)
371              ISW = NEWIDX(JPSW,NEXT)
372              ISE = NEWIDX(JPSE,NEXT)
373C
374C             Test if any of the four neighbouring points is missing.
375C
376                IF( (OLDWAVE(INW).EQ.PMISS) .OR.
377     X              (OLDWAVE(ISW).EQ.PMISS) .OR.
378     X              (OLDWAVE(INE).EQ.PMISS) .OR.
379     X              (OLDWAVE(ISE).EQ.PMISS) ) THEN
380                ENDIF
381
382              IF( (INW.EQ.0) .OR. (OLDWAVE(INW).EQ.PMISS) .OR.
383     X            (ISW.EQ.0) .OR. (OLDWAVE(ISW).EQ.PMISS) .OR.
384     X            (INE.EQ.0) .OR. (OLDWAVE(INE).EQ.PMISS) .OR.
385     X            (ISE.EQ.0) .OR. (OLDWAVE(ISE).EQ.PMISS) ) THEN
386                IF( (OLDWAVE(INW).EQ.PMISS) .AND.
387     X              (OLDWAVE(ISW).EQ.PMISS) .AND.
388     X              (OLDWAVE(INE).EQ.PMISS) .AND.
389     X              (OLDWAVE(ISE).EQ.PMISS) ) THEN
390                ENDIF
391C
392C               If so, take nearest grid point value.
393C
394                DISNW = DISTNEW(JPDISNW,NEXT)
395                DISNE = DISTNEW(JPDISNE,NEXT)
396                DISSW = DISTNEW(JPDISSW,NEXT)
397                DISSE = DISTNEW(JPDISSE,NEXT)
398C
399                IF( (DISNW.LE.DISNE).AND.
400     X              (DISNW.LE.DISSW).AND.
401     X              (DISNW.LE.DISSE)) THEN
402                  INDEX = INW
403C
404                ELSE IF( (DISNE.LE.DISSW).AND.
405     X                   (DISNE.LE.DISSE) ) THEN
406                  INDEX = INE
407C
408                ELSE IF( (DISSW.LE.DISSE) ) THEN
409                  INDEX = ISW
410C
411                ELSE
412                  INDEX = ISE
413                ENDIF
414C
415                IF(INDEX.EQ.0.OR.(OLDWAVE(INDEX).EQ.PMISS)) THEN
416C
417C                 Nearest point is missing
418C
419                  NEWWAVE(NEXTWV) = PMISS
420C
421                ELSE
422                  NEWWAVE(NEXTWV) =  OLDWAVE(INDEX)
423                ENDIF
424C
425              ELSE
426C
427C               Use bi-linear interpolation from four
428C               neighbouring sea points (possible conv REAL*4 to REAL*8)
429C
430                DI1N = DISTNEW(JPNW,NEXT)
431                DI2N = DISTNEW(JPNE,NEXT)
432                DI1S = DISTNEW(JPSW,NEXT)
433                DI2S = DISTNEW(JPSE,NEXT)
434                DK1  = DISTNEW(JPN,NEXT)
435                DK2  = DISTNEW(JPS,NEXT)
436C
437                NW_PT = OLDWAVE(INW)
438                NE_PT = OLDWAVE(INE)
439                SW_PT = OLDWAVE(ISW)
440                SE_PT = OLDWAVE(ISE)
441                IF( .NOT. LDIREC ) THEN
442                  U1 = NW_PT*DI2N + NE_PT*DI1N
443                  U2 = SW_PT*DI2S + SE_PT*DI1S
444                  NEWWAVE(NEXTWV) = U1*DK2 + U2*DK1
445                ELSE
446C
447C                 Fields representing a 'direction': resolve into
448C                 components and interpolate each separately.
449C
450#if (defined REAL_8)
451#define __ATAN2 DATAN2
452#define __COS   DCOS
453#define __SIN   DSIN
454#else
455#define __ATAN2 ATAN2
456#define __COS   COS
457#define __SIN   SIN
458#endif
459                  CNW_PT = __COS(NW_PT*RAD)
460                  CNE_PT = __COS(NE_PT*RAD)
461                  CSW_PT = __COS(SW_PT*RAD)
462                  CSE_PT = __COS(SE_PT*RAD)
463                  SNW_PT = __SIN(NW_PT*RAD)
464                  SNE_PT = __SIN(NE_PT*RAD)
465                  SSW_PT = __SIN(SW_PT*RAD)
466                  SSE_PT = __SIN(SE_PT*RAD)
467                  C1 = CNW_PT*DI2N + CNE_PT*DI1N
468                  C2 = CSW_PT*DI2S + CSE_PT*DI1S
469                  CC = C1*DK2 + C2*DK1
470                  S1 = SNW_PT*DI2N + SNE_PT*DI1N
471                  S2 = SSW_PT*DI2S + SSE_PT*DI1S
472                  SS = S1*DK2 + S2*DK1
473                  IF( SS.LT.0.0 ) THEN
474                    NEWWAVE(NEXTWV) = __ATAN2(SS,CC)/RAD + 360.0
475                  ELSE
476                    NEWWAVE(NEXTWV) = __ATAN2(SS,CC)/RAD
477                  ENDIF
478#undef __SIN
479#undef __COS
480#undef __ATAN2
481                ENDIF
482              ENDIF
483            ENDIF
484          ENDDO
485C
486        ENDIF
487        ISTART = ISTART + NUMNEW(NROW)
488      ENDDO
489C
490C     -----------------------------------------------------------------|
491C*    Section 9. Closedown.
492C     -----------------------------------------------------------------|
493C
494  900 CONTINUE
495      RETURN
496      END
497
498