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 WVQLIN2(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE,
13     X                         NEWWAVE,NORTH,WEST,KPARAM,PMISS)
14C
15C---->
16C*****WVQLIN2*
17C
18C     PURPOSE
19C     -------
20C
21C     To interpolate a wave field quasi-regular latitude
22C     longitude grid to a regular latitude/longitude grid.
23C
24C
25C     INTERFACE
26C     ---------
27C
28C     IRET = WVQLIN2(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE,
29C    X               NEWWAVE,NORTH,WEST,KPARAM,PMISS)
30C
31C     Input arguments
32C     ---------------
33C
34C     KNUM    - No. of meridians from North to South pole (input field)
35C     NUMPTS  - Array giving number of points along each latitude
36C               (empty latitudes have entry 0)
37C     KE_W    - First dimension of new array
38C               = Number of points E-W in new grid
39C     KN_S    - Second dimension of new array
40C               = Number of points N-S in new grid
41C     RESON   - Output grid resolution (degrees)
42C     OLDWAVE - Original wave field
43C     NORTH   - Input and output grid northernmost latitude (degree)
44C     WEST    - Input and output grid westernmost  longitude (degree)
45C     KPARAM  - Field parameter code
46C     PMISS   - Missing value indicator
47C
48C     Output arguments
49C     ----------------
50C
51C     NEWWAVE - New wave field
52C
53C     Function returns 0 if the interpolation was OK.
54C
55C
56C     METHOD
57C     ------
58C
59C     A bi-linear interpolation is only done if all four neighbouring
60C     points are sea.
61C     The nearest neighbouring grid point value is taken if one of the
62C     neighbouring points is missing.
63C     The interpolated point value is set to the 'missing' indicator
64C     value if the point is outside the input grid.
65C
66C     A field representing a 'direction' is resolved into components
67C     each of which is interpolated separately.
68C
69C     The output area is the same as the input area.
70C
71C
72C     EXTERNALS
73C     ---------
74C
75C     WVQLID2 - Calculate the indices of neighbouring points
76C     INTLOG  - Log error message.
77C
78C
79C     REFERENCE
80C     ---------
81C
82C     Based on:
83C         SUBROUTINE INTERPOLATE
84C         Peter Janssen     ECMWF    September 1995
85C     and:
86C         SUBROUTINE EXPOINT
87C         Heinz Gunther     ECMWF    December  1989
88C
89C
90C     Author.
91C     -------
92C
93C     J.D.Chambers      ECMWF    September 1996
94C
95C     J.D.Chambers      ECMWF    September 1998
96C     Modified to handle subarea input/output.
97C     The output area is the same as the input area.
98C
99C     J.D.Chambers      ECMWF    November 2000
100C     Modified to input grid resolution to be different from the
101C     output grid resolution.
102C
103C----<
104C
105      IMPLICIT NONE
106C
107#include "parim.h"
108C
109C     Parameters
110C
111      INTEGER JPROUTINE, JPLLMAX, JPNMOUT,
112     X        JPNW, JPNE, JPSW, JPSE, JPN, JPS,
113     X        JPDISNW, JPDISNE, JPDISSW, JPDISSE
114      PARAMETER( JPROUTINE = 19400 )
115      PARAMETER( JPLLMAX   =  1801 ) ! allow up to 0.1 degree resolution
116      PARAMETER( JPNMOUT   =  1800 ) ! allow up to 0.1 degree resolution
117      PARAMETER( JPNW=1, JPNE=2, JPSW=3, JPSE=4, JPN=5, JPS=6,
118     X           JPDISNW=7, JPDISNE=8, JPDISSW=9, JPDISSE=10 )
119C
120C     Arguments
121C
122      INTEGER KNUM, NUMPTS(*), KE_W, KN_S, KPARAM
123      REAL RESON, NORTH, WEST, PMISS
124      REAL RNS !FIXME: Unused! Difference in degrees in NS disrection
125      REAL OLDWAVE(*), NEWWAVE(*)
126C
127C     Local variables
128C
129      INTEGER NLAT, K, I, NEWROW, NEWCOL, I1N, I1S, I_N, I_S, I2N, I2S
130      REAL*4 DELONGN, DELONGS, DELAT, DI1N, DI1S, DI2N, DI2S
131      REAL*4 DIST_N, DIST_S
132      REAL*4 DISNW, DISNE, DISSW, DISSE
133      REAL*4 U1, U2, ZXIN, ZXIS
134      REAL LON(JPNMOUT*2), LAT(JPNMOUT+1)
135      LOGICAL INGRID, LDIREC
136      INTEGER NW_PT, NE_PT, SW_PT, SE_PT
137      REAL*4 NWVALUE, NEVALUE, SWVALUE, SEVALUE
138      REAL*4 RAD
139      REAL*4 CNW_PT, CNE_PT, CSW_PT, CSE_PT
140      REAL*4 SNW_PT, SNE_PT, SSW_PT, SSE_PT
141      REAL*4 C1, C2, S1, S2, CC, SS
142C
143      DATA RAD/0.017453293/
144C
145      REAL OLDLATS(JPLLMAX), OLDEAST
146      REAL OLDNS
147      INTEGER ISTART, IEND
148      INTEGER LOCATE
149      INTEGER NEXT, NUMNEW(JPLLMAX)
150      INTEGER IRET
151      INTEGER XKNUM,XKE_W,XKN_S
152      DATA XKNUM/-1/,XKE_W/-1/,XKN_S/-1/
153      REAL XRESON,XNORTH,XWEST
154      DATA XRESON/-999.0/,XNORTH/-999.0/,XWEST/-999.0/
155      INTEGER IFIRST, NEWIDX(4,3600*1801)
156      DATA IFIRST/0/
157      REAL*4 DISTNEW(10,3600*1801)
158C
159      INTEGER IOS
160      INTEGER OLDKNUM
161      DATA OLDKNUM/0/
162      INTEGER INDEX(JPLLMAX)
163      REAL DEPS
164      DATA DEPS/0.00005/
165C     Inline function - tests TRUE if A is a missing value
166C     EPS      Tolerance on missing data flag
167      REAL EPS
168      DATA EPS/0.1/
169      LOGICAL LSW,LSE,LNE,LNW,LTEMP
170C
171      SAVE IFIRST,XKNUM,XKE_W,XKN_S,XRESON,XNORTH,XWEST,NEWIDX,DISTNEW
172C
173C     Externals
174C
175      INTEGER WVQLID2
176      LOGICAL IS_WAVE_DIRECTION
177      EXTERNAL WVQLID2, IS_WAVE_DIRECTION
178C
179C     -----------------------------------------------------------------|
180C*    Section 1. Initalisation.
181C     -----------------------------------------------------------------|
182C
183  100 CONTINUE
184C
185      WVQLIN2 = 0
186      CALL INTLOG(JP_DEBUG,
187     X   'WVQLIN2: Wave interpolation requested.',JPQUIET)
188C
189C     Check reduced latitude/longitude grid specification
190C
191      IF( KNUM.GT.JPLLMAX ) THEN
192        CALL INTLOG(JP_ERROR,
193     X    'WVQLIN2: Number of latitudes in input lat/long grid = ',KNUM)
194        CALL INTLOG(JP_ERROR,
195     X    'WVQLIN2: And is greater than allowed maximum = ',JPLLMAX)
196        WVQLIN2 = JPROUTINE + 1
197        GOTO 900
198      ENDIF
199C
200C     Ensure working array dimensions are adequate for required output.
201C
202      IF( (KE_W.GT.JPNMOUT*2).OR.(KN_S.GT.JPNMOUT+1) ) THEN
203        CALL INTLOG(JP_ERROR,
204     X    'WVQLIN2: Internal array dimensions are too small',JPQUIET)
205        CALL INTLOG(JP_ERROR,
206     X    'WVQLIN2: for given lat/long output field.',JPQUIET)
207        WVQLIN2 = JPROUTINE + 3
208        GOTO 900
209      ENDIF
210C
211C     Set up index for latitude lines in the input reduced lat/lon array
212C
213      INDEX(1) = 0
214      DO K = 2, KNUM
215        INDEX(K) = INDEX(K-1) + NUMPTS(K-1)
216      ENDDO
217
218      DO K = 1, KN_S
219        LAT(K)     = NORTH - FLOAT(K-1)*RESON
220        NUMNEW(K)  = KE_W
221      ENDDO
222C
223C     Calculate latitudes and longitudes of output grid points.
224C
225      DO I = 1,KE_W
226        LON(I) = FLOAT(I-1)*RESON
227      ENDDO
228C
229C     Calculate latitudes for input (old) grid.
230C     (Input and output area are same)
231C
232      OLDNS =  RNS/REAL(KNUM-1)
233      DO K = 1, KNUM
234        OLDLATS(K) = NORTH - FLOAT(K-1)*OLDNS
235      ENDDO
236C
237C     Find first and last latitudes in input grid which are not empty
238C
239      ISTART = 1
240      DO K = 1, KNUM/2
241        IF( NUMPTS(K).EQ.0 ) ISTART = ISTART + 1
242      ENDDO
243C
244      IEND = KNUM
245      DO K = KNUM, (KNUM/2) + 1, -1
246        IF( NUMPTS(K).EQ.0 ) IEND = IEND - 1
247      ENDDO
248C
249C     For the regular output grid, calculate its latitudes and fill
250C     in the number of points along each latitude.
251C
252C
253C     Setup East depending on whether or not there is wrap-around
254C
255      IF( ABS((KE_W*RESON)-360.0).LT.DEPS ) THEN
256        OLDEAST = WEST + KE_W*RESON
257      ELSE
258        OLDEAST = WEST + (KE_W-1)*RESON
259      ENDIF
260C
261C     Initialise all points with 'missing data' indicator
262C
263      NEXT  = 0
264      DO NEWROW = 1, KN_S
265        DO NEWCOL = 1, KE_W
266          NEXT = NEXT + 1
267          NEWWAVE(NEXT) = PMISS
268        ENDDO
269      ENDDO
270C
271C     Wave direction parameters need special handling
272      LDIREC = IS_WAVE_DIRECTION(KPARAM)
273C
274C     -----------------------------------------------------------------|
275C*    Section 2. Interpolation.
276C     -----------------------------------------------------------------|
277C
278  200 CONTINUE
279C
280C     Only calculate the indices on the first time through
281C
282      IF( (IFIRST.EQ.0 )     .OR.
283     X    ( KNUM.NE.XKNUM )  .OR.
284     X    ( KE_W.NE.XKE_W )  .OR.
285     X    ( KN_S.NE.XKN_S )  .OR.
286     X    ( RESON.NE.XRESON ).OR.
287     X    ( NORTH.NE.XNORTH ).OR.
288     X    ( WEST.NE.XWEST ) ) THEN
289C
290        IRET =  WVQLID2(KNUM,NUMPTS,OLDLATS,WEST,OLDEAST,
291     X                  KN_S,NUMNEW,LAT,NEWIDX,DISTNEW)
292        IFIRST = 1
293        XKNUM  = KNUM
294        XKE_W  = KE_W
295        XKN_S  = KN_S
296        XRESON = RESON
297        XNORTH = NORTH
298        XWEST  = WEST
299      ENDIF
300C
301      DELAT = OLDNS
302      NEXT  = 0
303C
304      DO 220 NEWROW = 1, KN_S
305C
306C       Set up the distance between new row and the old rows to N and S.
307C
308        NLAT  = NINT((NORTH - LAT(NEWROW))/DELAT) + 1
309        I_N   = NLAT
310        I_S   = MIN(I_N+1,KNUM*2)
311        DIST_N  = ((NORTH - DELAT*(I_N-1)) - LAT(NEWROW)) / DELAT
312        DIST_S  = 1.0 - DIST_N
313C
314C       Check if the new interpolated row lies between 2 old rows which
315C       have data points
316C
317        IF( (NUMPTS(I_N).GT.0).AND.(NUMPTS(I_S).GT.0) ) THEN
318C
319C         Yes, use the calculated indices
320C
321          DO 210 NEWCOL = 1, KE_W
322C
323            NEXT = NEXT + 1
324C
325            NW_PT = NEWIDX(JPNW,NEXT)
326            NE_PT = NEWIDX(JPNE,NEXT)
327            SW_PT = NEWIDX(JPSW,NEXT)
328            SE_PT = NEWIDX(JPSE,NEXT)
329            NWVALUE = OLDWAVE(NW_PT)
330            NEVALUE = OLDWAVE(NE_PT)
331            SWVALUE = OLDWAVE(SW_PT)
332            SEVALUE = OLDWAVE(SE_PT)
333C
334C           Test if any one of the four neighbouring points is missing.
335C
336              LTEMP = .FALSE.
337              LSE = .FALSE.
338             IF(ABS(SEVALUE-PMISS) .LT. EPS)   LSE = .TRUE.
339              LSW = .FALSE.
340             IF(ABS(SWVALUE-PMISS) .LT. EPS)   LSW = .TRUE.
341              LNE = .FALSE.
342             IF(ABS(NEVALUE-PMISS) .LT. EPS)   LNE = .TRUE.
343              LNW = .FALSE.
344             IF(ABS(NWVALUE-PMISS) .LT. EPS)   LNW = .TRUE.
345C
346            IF( (NW_PT.EQ.0 .OR. LNW) .OR.
347     X          (NE_PT.EQ.0 .OR. LNE) .OR.
348     X          (SW_PT.EQ.0 .OR. LSW) .OR.
349     X          (SE_PT.EQ.0 .OR. LSE) ) THEN
350C             If so, take nearest grid point value.
351C
352              DISNW = DISTNEW(JPDISNW,NEXT)
353              DISNE = DISTNEW(JPDISNE,NEXT)
354              DISSW = DISTNEW(JPDISSW,NEXT)
355              DISSE = DISTNEW(JPDISSE,NEXT)
356C
357              IF( (DISNW.LE.DISNE).AND.
358     X            (DISNW.LE.DISSW).AND.
359     X            (DISNW.LE.DISSE) ) THEN
360cs                NEWWAVE(NEXT) = NWVALUE
361                LOCATE = NW_PT
362                LTEMP  = LNW
363C
364              ELSE IF( (DISNE.LE.DISSW).AND.
365     X                 (DISNE.LE.DISSE) ) THEN
366cs                NEWWAVE(NEXT) = NEVALUE
367                LOCATE = NE_PT
368                LTEMP  = LNE
369C
370              ELSE IF( (DISSW.LE.DISSE) ) THEN
371cs                NEWWAVE(NEXT) = SWVALUE
372                LOCATE = SW_PT
373                LTEMP  = LSW
374C
375              ELSE
376cs                NEWWAVE(NEXT) = SEVALUE
377                LOCATE = SE_PT
378                LTEMP  = LSE
379C
380              ENDIF
381
382              IF( .NOT.LTEMP.AND.LOCATE.NE.0) THEN
383c              IF( LOCATE.NE.0) THEN
384                  NEWWAVE(NEXT) = OLDWAVE(LOCATE)
385              ENDIF
386
387C
388            ELSE
389C
390C             Use bi-linear interpolation from four
391C             neighbouring sea points.
392C
393              DI1N = DISTNEW(JPNW,NEXT)
394              DI2N = DISTNEW(JPNE,NEXT)
395              DI1S = DISTNEW(JPSW,NEXT)
396              DI2S = DISTNEW(JPSE,NEXT)
397              DIST_N  = DISTNEW(JPN,NEXT)
398              DIST_S  = DISTNEW(JPS,NEXT)
399C
400              NW_PT = OLDWAVE(NW_PT)
401              NE_PT = OLDWAVE(NE_PT)
402              SW_PT = OLDWAVE(SW_PT)
403              SE_PT = OLDWAVE(SE_PT)
404C
405              IF( .NOT. LDIREC ) THEN
406                U1 = NWVALUE*DI2N + NEVALUE*DI1N
407                U2 = SWVALUE*DI2S + SEVALUE*DI1S
408                NEWWAVE(NEXT) = U1*DIST_S + U2*DIST_N
409C
410              ELSE
411C
412C               Fields representing a 'direction': resolve into
413C               components and interpolate each separately.
414C
415                CNW_PT = COS(NWVALUE*RAD)
416                CNE_PT = COS(NEVALUE*RAD)
417                CSW_PT = COS(SWVALUE*RAD)
418                CSE_PT = COS(SEVALUE*RAD)
419                SNW_PT = SIN(NWVALUE*RAD)
420                SNE_PT = SIN(NEVALUE*RAD)
421                SSW_PT = SIN(SWVALUE*RAD)
422                SSE_PT = SIN(SEVALUE*RAD)
423                C1 = CNW_PT*DI2N + CNE_PT*DI1N
424                C2 = CSW_PT*DI2S + CSE_PT*DI1S
425                CC = C1*DIST_S + C2*DIST_N
426                S1 = SNW_PT*DI2N + SNE_PT*DI1N
427                S2 = SSW_PT*DI2S + SSE_PT*DI1S
428                SS = S1*DIST_S + S2*DIST_N
429                IF( SS.LT.0.0 ) THEN
430                  NEWWAVE(NEXT) = ATAN2(SS,CC)/RAD + 360.0
431                ELSE
432                  NEWWAVE(NEXT) = ATAN2(SS,CC)/RAD
433                ENDIF
434              ENDIF
435            ENDIF
436  210     CONTINUE
437        ELSE
438           NEXT = NEXT + KE_W
439        ENDIF
440C
441  220 CONTINUE
442C
443C     -----------------------------------------------------------------|
444C*    Section 9. Closedown.
445C     -----------------------------------------------------------------|
446C
447  900 CONTINUE
448C
449      IF( WVQLIN2.EQ.0 ) THEN
450        CALL INTLOG(JP_DEBUG,
451     X    'WVQLIN2: Wave interpolation successful.',JPQUIET)
452      ELSE
453        CALL INTLOG(JP_DEBUG,
454     X    'WVQLIN2: Wave interpolation failed.',JPQUIET)
455      ENDIF
456c     call abort()
457C
458      RETURN
459      END
460