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 WV2DINT(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE,
12     X                         NEWWAVE,NORTH,WEST,KNSPEC,PMISS,RNS)
13C
14C---->
15C*****WV2DINT*
16C
17C     PURPOSE
18C     -------
19C
20C     To interpolate a wave field quasi-regular latitude
21C     longitude grid to a regular latitude/longitude grid.
22C
23C
24C     INTERFACE
25C     ---------
26C
27C     IRET = WV2DINT(KNUM,NUMPTS,KE_W,KN_S,RESON,OLDWAVE,
28C                    NEWWAVE,NORTH,WEST,KNSPEC,PMISS,RNS)
29C
30C     Input arguments
31C     ---------------
32C
33C     KNUM    - No. of meridians from North to South pole (input field)
34C     NUMPTS  - Array giving number of points along each latitude
35C               (empty latitudes have entry 0)
36C     KE_W    - First dimension of new array
37C               = Number of points E-W in new grid
38C     KN_S    - Second dimension of new array
39C               = Number of points N-S in new grid
40C     RESON   - Output grid resolution (degrees)
41C     OLDWAVE - Original wave field
42C     NORTH   - Output grid northernmost latitude (degree)
43C     WEST    - Output grid westernmost  longitude (degree)
44C     KNSPEC  - Number of 2D spectra values at each wave point
45C     PMISS   - Missing data value
46C     RNS     - Difference in degrees in NS disrection
47C
48C     Output arguments
49C     ----------------
50
51C     NEWWAVE - New wave field
52C
53C     Function returns 0 if the interpolation was OK.
54C
55C
56C     METHOD
57C     ------
58C
59C     The nearest neighbouring grid point value is assigned to
60C     the interpolated point.
61C
62C     The output (sub) is the same as the input area.
63C
64C
65C     EXTERNALS
66C     ---------
67C
68C     INTLOG  - Log error message.
69C
70C
71C     REFERENCE
72C     ---------
73C
74C     Based on:
75C         SUBROUTINE INTERPOLATE
76C         Peter Janssen     ECMWF    September 1995
77C     and:
78C         SUBROUTINE EXPOINT
79C         Heinz Gunther     ECMWF    December  1989
80C
81C
82C     Author.
83C     -------
84C
85C     J.D.Chambers      ECMWF    November 1996
86C
87C     J.D.Chambers      ECMWF    September 1998
88C     Modified to handle subarea input/output.
89C     The output (sub) is the same as the input area.
90C
91C
92C----<
93C
94      IMPLICIT NONE
95C
96#include "parim.h"
97C
98C     Parameters
99      INTEGER JPROUTINE
100      PARAMETER ( JPROUTINE = 19410 )
101      INTEGER JPLLMAX
102      PARAMETER ( JPLLMAX = 1801 )
103C                             `--> allow upto 0.1 degree resolution for
104C                                  Mediterranean, 0.1 for Global
105      INTEGER JPNMOUT
106      PARAMETER ( JPNMOUT = 1800 )
107C                             `--> allow upto 0.1 degree resolution for
108C                                  Mediterranean, 0.1 for Global
109C
110C     Subroutine arguments
111C
112      INTEGER KNUM, NUMPTS, KE_W, KN_S, KNSPEC
113      DIMENSION NUMPTS(*)
114      REAL RESON, OLDWAVE, NEWWAVE, NORTH, WEST, PMISS
115      REAL RNS
116      DIMENSION OLDWAVE(*)
117      DIMENSION NEWWAVE(KE_W,KN_S)
118C
119C     Local arguments
120C
121      INTEGER NLAT, K, NEWCOL, NEWROW
122      INTEGER I_N, I_S
123      REAL DELAT
124      REAL DIST_N, DIST_S
125      REAL LAT
126      DIMENSION LAT(JPNMOUT+1)
127      INTEGER INDEX
128      DIMENSION INDEX(JPLLMAX)
129C
130      REAL OLDLATS(JPLLMAX), OLDEAST
131      INTEGER NEXT, NUMNEW(JPLLMAX)
132      INTEGER IRET, W251IDX
133      EXTERNAL W251IDX
134      INTEGER XKNUM,XKE_W,XKN_S
135      DATA XKNUM/-1/,XKE_W/-1/,XKN_S/-1/
136      REAL XRESON,XNORTH,XWEST
137      DATA XRESON/-999.0/,XNORTH/-999.0/,XWEST/-999.0/
138      INTEGER IFIRST, NEWIDX(3600*1801)
139      DATA IFIRST/0/
140      REAL DEPS
141      DATA DEPS/0.00005/
142C
143      SAVE IFIRST,XKNUM,XKE_W,XKN_S,XRESON,XNORTH,XWEST,NEWIDX
144C
145C
146C ---------------------------------------------------------------------
147C*    Section 1. Initalisation.
148C ---------------------------------------------------------------------
149C
150  100 CONTINUE
151C
152      WV2DINT = 0
153      CALL INTLOG(JP_DEBUG,
154     X   'WV2DINT: Wave interpolation requested.',JPQUIET)
155C
156C     Check only new-style 2D wave spectra (parameter 251)
157C
158      IF( KNSPEC.GT.1 ) THEN
159        CALL INTLOG(JP_ERROR,
160     X    'WV2DINT: Only single-value 2D spectra field allowed',JPQUIET)
161        CALL INTLOG(JP_ERROR,
162     X    'WV2DINT: Value given = ',KNSPEC)
163        WV2DINT = JPROUTINE + 1
164        GOTO 900
165      ENDIF
166C
167C     Check reduced latitude/longitude grid specification
168C
169      IF( KNUM.GT.JPLLMAX ) THEN
170        CALL INTLOG(JP_ERROR,
171     X    'WV2DINT: Number of latitudes in input lat/long grid = ',KNUM)
172        CALL INTLOG(JP_ERROR,
173     X    'WV2DINT: And is greater than allowed maximum = ',JPLLMAX)
174        WV2DINT = JPROUTINE + 2
175        GOTO 900
176      ENDIF
177C
178C     Ensure working array dimensions are adequate for required output.
179C
180      IF( (KE_W.GT.JPNMOUT*2).OR.(KN_S.GT.JPNMOUT+1) ) THEN
181        CALL INTLOG(JP_ERROR,
182     X    'WV2DINT: Internal array dimensions are too small',JPQUIET)
183        CALL INTLOG(JP_ERROR,
184     X    'WV2DINT: for given lat/long output field.',JPQUIET)
185        WV2DINT = JPROUTINE + 3
186        GOTO 900
187      ENDIF
188C
189C     Set up index for latitude lines in the input reduced lat/lon array
190C
191      INDEX(1) = 0
192      DO K = 2, KNUM
193        INDEX(K) = INDEX(K-1) + NUMPTS(K-1)
194      ENDDO
195C
196C     Calculate latitudes and longitudes of output grid points.
197C
198C
199      DO K = 1, KN_S
200        LAT(K) = NORTH - FLOAT(K-1)*RESON
201        NUMNEW(K)  = KE_W ! 0
202      ENDDO
203
204      DELAT  = RNS/(KNUM-1)
205
206      DO k=1,KNUM
207          OLDLATS(K) = NORTH - FLOAT(K-1)*DELAT
208      END DO
209C
210      IF( ABS((KE_W*RESON)-360.0).LT.DEPS ) THEN
211        OLDEAST = WEST + KE_W*RESON
212      ELSE
213        OLDEAST = WEST + (KE_W-1)*RESON
214      ENDIF
215C
216C     Initialise all points with 'missing data' indicator
217C
218      DO NEWROW = 1, KN_S
219        DO NEWCOL = 1, KE_W
220          NEWWAVE(NEWCOL,NEWROW) = PMISS
221        ENDDO
222      ENDDO
223C
224C ---------------------------------------------------------------------
225C*    Section 2. Interpolation.
226C ---------------------------------------------------------------------
227C
228  200 CONTINUE
229C
230C     Only calculate the indices on the first time through
231C
232      IF( (IFIRST.EQ.0 )     .OR.
233     X    ( KNUM.NE.XKNUM )  .OR.
234     X    ( KE_W.NE.XKE_W )  .OR.
235     X    ( KN_S.NE.XKN_S )  .OR.
236     X    ( RESON.NE.XRESON ).OR.
237     X    ( NORTH.NE.XNORTH ).OR.
238     X    ( WEST.NE.XWEST ) ) THEN
239        IRET = W251IDX(KNUM,NUMPTS,OLDLATS,WEST,OLDEAST,
240     X                 KN_S,NUMNEW,LAT,NEWIDX)
241        IFIRST = 1
242        XKNUM  = KNUM
243        XKE_W  = KE_W
244        XKN_S  = KN_S
245        XRESON = RESON
246        XNORTH = NORTH
247        XWEST  = WEST
248      ENDIF
249C
250C      DELAT  = 180.0/(KNUM-1)
251      DELAT  = RNS/(KNUM-1)
252      NEXT = 0
253C
254      DO 220 NEWROW = 1, KN_S
255C
256C       Set up the distance between new row and the old rows to N and S.
257C
258        NLAT   = NINT((NORTH - LAT(NEWROW))/DELAT) + 1
259        I_N    = NLAT
260        I_S    = MIN(I_N+1,KNUM*2)
261        DIST_N = ((NORTH - DELAT*(I_N-1)) - LAT(NEWROW)) / DELAT
262        DIST_S = 1.0 - DIST_N
263C
264C       Check if the new interpolated row lies between 2 old rows which
265C       have data points
266C
267        IF( (NUMPTS(I_N).GT.0).AND.(NUMPTS(I_S).GT.0) ) THEN
268C
269C         Yes, use the calculated indices
270C
271          DO 210 NEWCOL = 1, KE_W
272C
273            NEXT = NEXT + 1
274C
275            NEWWAVE(NEWCOL,NEWROW) = OLDWAVE(NEWIDX(NEXT))
276C
277  210     CONTINUE
278        Else
279           NEXT = NEXT + KE_W
280        ENDIF
281C
282  220 CONTINUE
283C
284C ---------------------------------------------------------------------
285C*    Section 9. Closedown.
286C ---------------------------------------------------------------------
287C
288  900 CONTINUE
289C
290      IF( WV2DINT.EQ.0 ) THEN
291        CALL INTLOG(JP_DEBUG,
292     X    'WV2DINT: Wave interpolation successful.',JPQUIET)
293      ELSE
294        CALL INTLOG(JP_DEBUG,
295     X    'WV2DINT: Wave interpolation failed.',JPQUIET)
296      ENDIF
297C
298      RETURN
299      END
300