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 WV2DIDX(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST,
12     X                         KNEWNUM,NUMNEW,NEWLATS,NEWWAVE)
13C
14C---->
15C*****WV2DIDX*
16C
17C     PURPOSE
18C     -------
19C
20C     Determines which nearest-neighbour values of an input global wave
21C     2D-spectra grid field to use for an output global wave 2D-spectra
22C     grid field.
23C
24C
25C     INTERFACE
26C     ---------
27C
28C     IRET = WV2DIDX(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST,
29C    X               KNEWNUM,NUMNEW,NEWLATS,NEWWAVE)
30C
31C     Input arguments
32C     ---------------
33C
34C     KOLDNUM - No. of meridians from North to South pole (input field)
35C     NUMOLD  - Array giving number of points along each latitude
36C               (empty latitudes have entry 0)
37C     OLDLATS - input field latitudes
38C     OLDWEST - western longitude of the input field (degrees)
39C     OLDEAST - eastern longitude of the input field (degrees)
40C
41C     KNEWNUM - No. of meridians from North to South pole (output field)
42C     NUMNEW  - Array giving number of points along each latitude
43C               (empty latitudes have entry 0)
44C     NEWLATS - output field latitudes
45C
46C     Output arguments
47C     ----------------
48C
49C     NEWWAVE - Indices of points to use
50C
51C     Function returns 0 if the interpolation was OK.
52C
53C
54C     METHOD
55C     ------
56C
57C     The index of the nearest neighbouring grid point value of
58C     the input field is assigned to each point of the output
59C     grid.
60C     The index is zero if the output grid point is 'missing'.
61C
62C     The input field can be regular or quasi-regular, and can be
63C     global or a subarea.
64C     The output field can be regular or quasi-regular.
65C
66C
67C     EXTERNALS
68C     ---------
69C
70C     INTLOG  - Log error message.
71C
72C
73C     REFERENCE
74C     ---------
75C
76C     None
77C
78C
79C     Author.
80C     -------
81C
82C     J.D.Chambers      ECMWF    October  1997
83C
84C
85C----<
86C
87      IMPLICIT NONE
88C
89#include "parim.h"
90C
91C     Parameters
92C
93      INTEGER JPROUTINE
94      PARAMETER ( JPROUTINE = 19410 )
95      INTEGER JPLLMAX
96      PARAMETER ( JPLLMAX = 1801 )
97C                             `--> allow upto 0.1 degree resolution
98      INTEGER JPNMOUT
99      PARAMETER ( JPNMOUT = 1800 )
100C                             `--> allow upto 0.1 degree resolution
101C
102C     Function arguments
103C
104      INTEGER KOLDNUM, NUMOLD, KNEWNUM, NUMNEW
105      DIMENSION NUMOLD(*), NUMNEW(*)
106      REAL OLDWEST, OLDEAST, OLDLATS, NEWLATS
107      DIMENSION OLDLATS(*), NEWLATS(*)
108      INTEGER NEWWAVE
109      DIMENSION NEWWAVE(*)
110C
111C     Local arguments
112C
113      INTEGER NEWCOL, NEWROW, INCOL, LOOP, NEXT, NUMMAX
114      INTEGER I_NW, I_SW, I_N, I_S, I_NE, I_SE
115      REAL*8 DELONGN, DELONGS, DIFF, RESOL, DINC, DEPS
116      REAL*8 DIST_NW, DIST_SW, DIST_NE, DIST_SE, DIST_N, DIST_S
117      REAL*4 DISNW, DISNE, DISSW, DISSE
118      REAL*8 ZXIN, ZXIS, LON, OWEST, OEAST
119      DIMENSION LON(JPNMOUT*2)
120      LOGICAL LINGNS, LINGWE
121      INTEGER INDEXI
122      DIMENSION INDEXI(JPLLMAX)
123      LOGICAL LGLOBAL, LFULLG, LGWRAP
124      DATA DEPS/0.00005/
125C
126C ---------------------------------------------------------------------
127C*    Section 1. Initalisation.
128C ---------------------------------------------------------------------
129C
130  100 CONTINUE
131C
132      WV2DIDX = 0
133      CALL INTLOG(JP_DEBUG,
134     X   'WV2DIDX: Wave interpolation requested.',JPQUIET)
135C
136C     Check latitude/longitude grid specification
137      IF( KOLDNUM.GT.JPLLMAX ) THEN
138        CALL INTLOG(JP_ERROR,
139     X    'WV2DIDX: Number of latitudes in input grid = ',KOLDNUM)
140        CALL INTLOG(JP_ERROR,
141     X    'WV2DIDX: And is greater than allowed maximum = ',JPLLMAX)
142        WV2DIDX = JPROUTINE + 1
143        GOTO 900
144      ENDIF
145C
146C     Set up INDEXI for latitude lines in input lat/lon array
147C
148
149      INDEXI(1) = 0
150      DO LOOP = 2, KOLDNUM
151        INDEXI(LOOP) = INDEXI(LOOP-1) + NUMOLD(LOOP-1)
152      ENDDO
153C
154      OEAST = OLDEAST
155      OWEST = OLDWEST
156
157C     Check whether the input field is global or a subarea
158C
159      NUMMAX = NUMOLD(1)
160      DO LOOP = 2, KOLDNUM
161        IF( NUMOLD(LOOP).GT.NUMMAX ) NUMMAX = NUMOLD(LOOP)
162      ENDDO
163      IF( NUMMAX.LE.0 ) THEN
164        CALL INTLOG(JP_ERROR,
165     X    'WV2DIDX: Input wave field has no points',JPQUIET)
166        WV2DIDX = JPROUTINE + 2
167        GOTO 900
168      ENDIF
169C
170      LGWRAP = ABS((OEAST-OWEST) - DBLE(360.0)).LT.DEPS
171      DINC = 360.0/DBLE(NUMMAX)
172      LFULLG = ABS((DMOD((OEAST-OWEST+DBLE(360.0)),DBLE(360.0))
173     X              +DINC)-DBLE(360.0)).LT.DEPS
174C
175      LGLOBAL = LFULLG.OR.LGWRAP
176C
177C ---------------------------------------------------------------------
178C*    Section 2. Interpolation.
179C ---------------------------------------------------------------------
180C
181  200 CONTINUE
182C
183      NEXT = 1
184C
185      DO 220 NEWROW = 1, KNEWNUM
186C
187C       Find old latitude rows to north and south of new latitude and
188C       set up the distance between new row and the old rows to N and S.
189C
190        DO LOOP = 1, KOLDNUM -1
191          IF( (NEWLATS(NEWROW).LE.OLDLATS(LOOP)) .AND.
192     X        (NEWLATS(NEWROW).GT.OLDLATS(LOOP+1)) ) THEN
193            I_N = LOOP
194            I_S = MIN(LOOP+1,KOLDNUM)
195            GOTO 205
196          ENDIF
197        ENDDO
198C
199        IF( NEWLATS(NEWROW).EQ.OLDLATS(KOLDNUM) ) THEN
200          I_N = KOLDNUM
201          I_S = KOLDNUM
202        ELSE
203          I_N = -1
204          I_S = -1
205        ENDIF
206C
207  205   CONTINUE
208        IF( I_N.NE.I_S) THEN
209          DIFF   = OLDLATS(I_N) - OLDLATS(I_S)
210          DIST_N = (OLDLATS(I_N) - NEWLATS(NEWROW)) / DIFF
211        ELSE
212          DIST_N = 1.0
213        ENDIF
214        DIST_S = 1.0 - DIST_N
215C
216C       Check if the new interpolated row lies between 2 old rows which
217C       have data points
218C
219        LINGNS = .FALSE.
220        IF( I_N.GT.0 ) THEN
221          IF( (NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).GT.0) ) THEN
222C
223C           Yes, so set up the grid increments.
224C
225C
226C           Note that the grid increments are different depending on
227C           whether or not the grid is global; global grids have an
228C           extra interval at the end (eg between 358.5 and 0.0), while
229C           subareas have outer limits defined by 'west' and 'east' and
230C           a given number of points in between.
231C
232            IF( LGLOBAL ) THEN
233              DELONGN = 360.0/DBLE(NUMOLD(I_N))
234              DELONGS = 360.0/DBLE(NUMOLD(I_S))
235            ELSE
236              DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/
237     X                 DBLE(NUMOLD(I_N))
238
239              DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
240     X                 +DINC) / DBLE(NUMOLD(I_N))
241              DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/
242     X                 DBLE(NUMOLD(I_S))
243              DELONGS = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
244     X                 +DINC) / DBLE(NUMOLD(I_S))
245            ENDIF
246c
247
248            LINGNS  = .TRUE.
249          ENDIF
250C
251C         The equator is given special treatment so that the northern
252C         and southern hemispheres for the parameter 250 can (later)
253C         be stitched together.
254C
255C         If the input field finishs at the equator and the output
256C         field has a line at the equator, use the input equator for
257C         interpolation.
258C
259          IF( NEWLATS(NEWROW).EQ.0.0 .AND.
260     x        ((NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).EQ.0)) ) THEN
261C
262C           Yes, so set up the grid increments.
263C
264            IF( LGLOBAL ) THEN
265              DELONGN = 360.0/DBLE(NUMOLD(I_N))
266            ELSE
267              DINC = DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))/
268     X                 DBLE(NUMOLD(I_N))
269              DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
270     X                 +DINC) / DBLE(NUMOLD(I_N))
271            ENDIF
272            DELONGS = DELONGN
273            DIST_N = 0.0
274            DIST_S = 1.0
275            I_S = I_N
276            LINGNS  = .TRUE.
277          ENDIF
278        ENDIF
279C
280C       Setup longitudes for the current row
281C
282        INCOL = NUMNEW(NEWROW)
283        IF( INCOL.NE.0 ) THEN
284             RESOL = 360.0/DBLE(INCOL)
285             DO LOOP = 1, INCOL
286               LON(LOOP) = DBLE(LOOP-1)*RESOL
287             ENDDO
288C
289          DO 210 NEWCOL = 1, INCOL
290C
291            LINGWE = (LON(NEWCOL).GT.OWEST).OR.(LON(NEWCOL).LT.OEAST)
292C
293            IF( LINGNS .AND. LINGWE ) THEN
294C
295C             If point to be interpolated is inside the input field ,
296C             work out distances from new point to four neighbours.
297C
298              ZXIN = DMOD( (LON(NEWCOL) + (360.0-OWEST) + DBLE(720.0)),
299     X                     DBLE(360.0))
300cs            ZXIN = LON(NEWCOL)
301              ZXIS = ZXIN
302              ZXIN = ZXIN/DELONGN + 1.0
303              ZXIS = ZXIS/DELONGS + 1.0
304              I_NW = INT(ZXIN)
305              I_NE = I_NW + 1
306              I_SW = INT(ZXIS)
307              I_SE = I_SW + 1
308C
309C             Allow wrap-around
310C
311C
312              I_NW = MOD(I_NW,NUMOLD(I_N))
313              I_NE = MOD(I_NE,NUMOLD(I_N))
314              I_SW = MOD(I_SW,NUMOLD(I_S))
315              I_SE = MOD(I_SE,NUMOLD(I_S))
316              IF( I_NW.EQ.0 ) I_NW = NUMOLD(I_N)
317              IF( I_SW.EQ.0 ) I_SW = NUMOLD(I_S)
318              IF( I_NE.EQ.0 ) I_NE = NUMOLD(I_N)
319              IF( I_SE.EQ.0 ) I_SE = NUMOLD(I_S)
320C
321C             Calculate distance from interpolated point to its neighbours
322C
323              DIST_NW = ZXIN - REAL(I_NW)
324              DIST_NE = 1.0 - DIST_NW
325              DIST_SW = ZXIS - REAL(I_SW)
326              DIST_SE = 1.0 - DIST_SW
327C
328C             Take nearest grid point value.
329C
330              DISNW = DIST_NW*DIST_NW + DIST_N*DIST_N
331              DISNE = DIST_NE*DIST_NE + DIST_N*DIST_N
332              DISSW = DIST_SW*DIST_SW + DIST_S*DIST_S
333              DISSE = DIST_SE*DIST_SE + DIST_S*DIST_S
334C
335              IF( (DISNW.LE.DISNE).AND.(DISNW.LE.DISSW).AND.
336     X            (DISNW.LE.DISSE) ) THEN
337C
338C               Use north-west neighbour
339C
340                NEWWAVE(NEXT) = INDEXI(I_N)+I_NW
341C
342              ELSE IF( (DISNE.LE.DISSW).AND.(DISNE.LE.DISSE) ) THEN
343C
344C               Use north-east neighbour
345C
346                NEWWAVE(NEXT) = INDEXI(I_N)+I_NE
347C
348              ELSE IF( DISSW.LE.DISSE ) THEN
349C
350C               Use south-west neighbour
351C
352                NEWWAVE(NEXT) = INDEXI(I_S)+I_SW
353              ELSE
354C
355C               Use south-east neighbour
356C
357                NEWWAVE(NEXT) = INDEXI(I_S)+I_SE
358              ENDIF
359C
360C           If the point to be interpolated is outside the input field,
361C           set it to 'land'.
362C
363            ELSE
364              NEWWAVE(NEXT) = 0
365            ENDIF
366C
367            NEXT = NEXT + 1
368C
369  210     CONTINUE
370C
371        ENDIF
372C
373  220 CONTINUE
374C
375C ---------------------------------------------------------------------
376C*    Section 9. Closedown.
377C ---------------------------------------------------------------------
378C
379  900 CONTINUE
380C
381      IF( WV2DIDX.EQ.0 ) THEN
382        CALL INTLOG(JP_DEBUG, 'WV2DIDX: successful.',JPQUIET)
383      ELSE
384        CALL INTLOG(JP_ERROR, 'WV2DIDX: failed.',JPQUIET)
385      ENDIF
386C
387      RETURN
388      END
389